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Hypersonic Nozzle Expansion of Air With 


Atom Recombination Present’ 


H. T. NAGAMATSU,* J. B. WORKMAN, ** anp R. E. SHEER, JR.*** 


General Electric Company 


Summary 


An experimental investigation on the expansion of high-tem- 
perature, high-pressure air to hypersonic flow Mach numbers in 
a conical nozzle of a hypersonic shock tunnel has been carried out. 
The equilibrium temperature and pressure ranges after the re- 
flected shock wave were 1400° to 6000°K and 100 to 1000 psia 
Static-pressure measurements, which are sensitive to the state 
of the gas, were made along the axis of the nozzle for different 
reservoir conditions. These results are compared with the cal- 
culated equilibrium and “‘frozen’’ data for the same geometry 
and initial reservoir conditions. 
® For reservoir pressures greater than 500 psia, the expansion of 
the air in the nozzle is essentially in equilibrium up to reservoir 
temperatures of about 4,500°K. For temperatures greater than 
4,000°K and pressures lower than 100 psia, the expanding air is 
almost frozen. At a given area ratio for the nozzle and reservoir 
pressure, the expansion process remains in equilibrium up to a 
certain reservoir temperature, and beyond this temperature the 
flow expansion deviates rapidly from the equilibrium process and 


approaches the frozen case. 


Introduction 


ITH THE GREAT INTEREST TODAY in the hyper- 

W sonic flight speeds, a number of test facilities 
have been constructed that attempt to simulate high- 
flight-Mach-number conditions. In most cases the 
high velocities are achieved by expanding heated air 
through a nozzle of large area ratio. If an attempt is 
made to simulate flight enthalpy in addition to Mach 
number, the expansion process must start from a res- 
ervoir of extremely high-temperature air. To achieve 
the velocity and enthalpy corresponding to Mach 20 
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at 100,000 ft altitude will require a reservoir tempera- 
ture of approximately 7,000°K for reasonable test- 
section pressures. 

For the expansion process from such a high tempera- 
ture the problems of atom recombination, chemical 
reactions, and vibrational relaxation will play dominant 
roles in the thermodynamics of the flow. Even at 
temperatures well below this relatively extreme case, 
a considerable amount of molecular dissociation will 
have taken place in the reservoir. For example, the 
oxygen in equilibrium air at 3,500°K and a pressure of 
one atmosphere is approximately 50 percent dissoci- 
ated.' Thus, atom recombination plays a large role 
in achieving the proper test-section conditions in hyper- 
sonic tunnel devices, and it is essential that an attempt 
be made to understand the expansion process. 

In general, to minimize the flow disturbances due to 
wall viscous effects for high flow Mach numbers, the 
hypersonic nozzle*: * attached to a shock tube is made as 
short as possible. This implies that the expansion 
process will be extremely rapid. If the residence time 
of the air in the nozzle is quite short, it becomes ques- 
tionable whether thermodynamic equilibrium can be 
maintained. While it can safely be assumed that the 
translational and rotational modes of energy storage 
will equilibrate in the time available, the modes of atom 
recombination, chemical reactions, and _ vibrational 
relaxation may be quite marginal in doing so. Since 
the lack of equilibrium of these modes can change the 
test-section pressure by as much as a factor of 10 from 
the equilibrium value, some determination of the re- 
action rates involved must be made. Therefore, in 
order to obtain significant data in a hypersonic test 
facility using an expansion nozzle, the exact state of the 
gas must be determined because of the possibility of 
nonequilibrium expansion. 

The results of investigations that were conducted in 
the nozzle and the test section of a hypersonic shock 
tunnel? at the General Electric Research Laboratory 
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Fic. 1. Isentropic pressure ratio as a function of geometric 
area ratio and reservoir temperature 


will be discussed in this paper. The departure of the 
expanded air from thermodynamic equilibrium was 
determined for a variety of reservoir pressures and tem- 
peratures. By this method, it was possible not only to 
calibrate the flow in the nozzle but also to gain some 
knowledge regarding the air recombination process. 

In recent years several theoretical papers!~'’ on the 
subject of dissociated air expanding in a nozzle have 
been prepared by assuming recombination rates. At 
th: present time there appears to be neither analytical 
nor experimental data adequate to describe the be- 
havior of air in the regime of interest. Therefore, an 
experimental investigation of a thermodynamic param- 
eter sensitive to departures from equilibrium was 
conducted. 

During the expansion process, the thermodynamic 
properties of air most affected by the deviation from 
equilibrium are the static pressure and temperature. 
Fortunately, the static pressure is the thermodynamic 
parameter that can be measured in a shock tunnel. 
It is possible to calculate'! the static pressure along an 
expansion nozzle for one-dimensional, isentropic, equi- 
librium flow and “‘frozen’’ flow.'? In this case, the 
““frozen”’ gas is defined to be a pseudo-equilibrium com- 
position wherein only translational and_ rotational 
energy modes equilibrate, all others remaining fixed 
at the reservoir condition. Presumably, if the nozzle 
flow is isentropic, these two calculations should bound 
the possible state of the test gas, and the measured 
static pressure will show the amount of departure from 
these two extremes. 

The above type of information is necessary for the 
calibration of a given hypersonic test facility. How- 
ever, in itself it is insufficient for determining the actual 
reaction rates in air because of the multitude of re- 
actions taking place simultaneously. Therefore, at 
present a correlation of the results with a chemical 
kinetic theory is difficult. Future study will include an 
investigation of the expansion of certain components of 
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air, such as nitrogen and oxygen, in the shock tunnel. 
It is hoped that eventually a complete description of 
the recombination and vibrational relaxation of air 
will be possible. The present results serve to show how 
a high stagnation temperature, hypersonic test facility 
may be calibrated, and will also give some qualitative 
information on dissociated air expansion processes in a 


hypersonic nozzle. 


Nozzle Analysis 


The method of calculating the flow in the equilibrium 
case for the nozzle is quite straightforward. The main 
flow is effectively one-dimensional and nearly free of 
viscous effects, as discussed in reference 2. This means 
that the thermodynamic state of the expanding gas can 
be determined as a function of the area ratio for a given 
reservoir condition 

The calculation then involves satisfying the conser- 
vation of mass and energy, and the equilibrium thermo- 
dynamic state relations. The necessary relationships 


are given by: 


plA p *A" (1) 
(n)12 +h Ino (2) 
p [(h,so) (:3) 


where p is the density, | is the velocity, .! is the cross 
sectional area, / is the enthalpy, and s is the entropy. 
The starred (*) quantities refer to conditions at the 
throat, while the zero subscript refers to initial reservoir 
conditions. In addition to the above equations, the 


throat condition of 
oe a* when A A* (4) 


must be satisfied, where a* is the frozen speed of sound" 


given by 
a g(h,so) (5 


> 


Of course, as implied by Eqs. (3) and (5), if and 
So are known, all other thermodynamic properties such 
as pressure and temperature are also known.'' The 
indicated functions / and g plus other similar relations 
among the thermodynamic properties are not analytic, 
but must be obtained from standard tabulated or 
plotted data," 1*-* 
tion of Feldman" was used above 2,000°K, while all 
other data were taken from Hirschfelder and Curtiss." 
The above relationships can be solved in an iterative 
While the 


present work was done by hand computation, the cal 


In the present work the compila 


manner to obtain the equilibrium solution. 


culations can be performed by an IBM 704 computer in 
a way devised by W. G. Browne, similar to that dis 
cussed in reference | 1. 

The calculation for a “‘frozen’’ expansion of air in 
the nozzle is actually simpler than that for an equi 
librium expansion. Since a frozen expansion!*:!’ has 
been defined as a gas change where no recombination or 
chemical reaction takes place and only translational 
and rotational energies are allowed to change during 
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HY PERSOXIC NOZZLE 


the flow process, any change in the enthalpy of the air 
is directly proportional to a change in the temperature. 
It must be borne in mind that by definition we are per 
mitting no recombination, vibrational relaxation, 
chemical change, etc., to occur. 

The assumption of the above flow model permits the 
calculation of an effective specific heat ratio defined by 


4 1 + (ZR/C,) (6) 


Ye 


where Z is the compressibility factor and C, is the spe- 
cific heat at constant volume. By definition, y, will 
be a constant. The above-defined y, has been cal- 
culated for air in the range of 1,500 to S,Q0O°K and 100 
The values of y, range from 1.4 to 1.67 
Once the 


to 3,000 psia. 
depending on the amount of gas dissociated. 
value for y, Eq. (6), is determined for the reservoir 
condition, the static pressure as a function of the area 
ratio can be calculated by use of the nozzle flow equa- 
tion for a perfect gas: 


& Me" [2/ (ve + 1)] (2/ Cre + 1)) 
yy 6 CS (P/Po)?/") [1 — (P/Po) 


\4 
(7) 

By the two methods outlined above, one for equi- 
librium and one for frozen expansion, it is possible to 
construct tabulations and graphs which show the static 
pressure as a function of area ratio, and reservoir pres- 
sure and temperature as indicated in Figs. 1-3. These 
calculated values should bound the experimental re- 
sults and aid in the interpretation of the process of 


expansion of the dissociated air. 
Experimental Technique 


The high-temperature, high-pressure reservoir air 
required for the expansion in the hypersonic nozzle for 
the investigation was generated by use of a combustion- 
driven shock tube.* In this tube it is possible to pro- 
duce, by a reflected shock wave from the entrance to 
the nozzle, compressed &nd heated air in equilibrium for 
a period of approximately | to 4 msec. Reservoir tem- 
peratures and pressures between 1,400°K, 1,000 psia 
and 6,000°K, 100 psia were investigated. 

The reservoir gas is expanded through a l-in. di- 
ameter throat with rounded edges of 1.5 1n. radius into a 
conical expansion nozzle with 30° included angle and 
an exit diameter of 24 in. 

Static pressure measurements, on the axis, were 
made in this nozzle over the diameter range of 7.5 in 
In this report the results at 7.5-in. and 1|2-in. 
These correspond 


to 24 in. 
diameters will be discussed in detail. 
to area ratios of 56.2 and 144 respectively. Mat- 
thews'* has reported on the effect that the geometric 
shape of the probe and the growth of the boundary 
layer along the probe will have on the measured static 
pressures. For the present probe this deviation will 
be a very small percentage when compared to the de- 
viation between frozen and equilibrium flow. These 
pressure measurements were made using a long slender 
probe with a diameter of '/s in. and a small-angle ogive 
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for the nose. The static pressure orifices were placed 
around the periphery of the probe 10.25 in. back from 
the tip. A barium titanate gage with a linear response 
of 12.5 mv./ psi was placed inside the probe to measure 
the static pressure. The rise time of this gage as in 
stalled in the probe was approximately 150 usec. The 
probe was attached to a 6-in. diameter sting some S ft 
long which was shock-isolated from the nozzle and the 
dump tank.” This prevented the mechanical vibra 
tion produced by the reflected shock wave at the nozzle 
entrance from being recorded by the probe 

In the larger-diameter portion of the nozzle a * ,-in 
diameter probe, with a gage having a 46 mv. psi out 
put and a 50) usec rise time, was also used. The static 
pressures varied from 0.08 to 0.60 psia with good agree 
ment between the two probes. Pressures as low as 
0.002 psia have been recorded with a 1-in. diameter 
probe and a pressure gage of 160 mv. /psia output 

A thin diaphragm was placed in the shock tube just 
ahead of the entrance to the nozzle in order that the 
nozzle and 200-ft.* dump tank could be evacuated to 
about 25 uw of mercury prior to the firing. This en 
abled the flow to start and to reach steady state in 
about 750 usec without any appreciable starting shock 
waves. 

Since water vapor is formed in the combustion of 
hydrogen and oxygen in the driver, the tube was care 
fully blown out and dried before each test. The reser 
voir gas is air of standard composition with a dew point 
of —S0°F. Thus, the gas is ‘dry’ with respect to 
equilibrium thermodynamic properties, but it may still 
have some effects on atom recombination and vibra 


tional relaxation. 


Results and Discussion 


A large number of tests were conducted in the shock 
tunnel with three different static pressure probes and 
an impact pressure probe which had a hemispherical 
nose shape. These experimental results are correlated 
with the calculated values for different reservoir pres 
sures, temperatures, and nozzle area ratios 

The impact pressure measurements along the axis 
of the nozzle are approximately the same whether the 
flow expansion process is equilibrium or frozen. These 
results agreed with the calculated values, as discussed 
in reference 12. Thus, the impact pressure measure 
ments, which are relatively easy to conduct, are not 
suitable for the determination of the state of the gas in 
the nozzle. The nozzle static pressure and the two 
dimensional wedge shock angles are calculated to be 
quite different for the two extremes of equilibrium and 
completely frozen flow expansions 

In Fig. 1 the axial static pressure results are presented 
of approximately 500 psia 
Also shown 


for a reservoir pressure /y 
and temperatures of 1,400° and 4,000°K 
are the calculated isentropic, one-dimensional flow 
values for the case of equilibrium and perfect gas (7 
1.4) expansion from the reservoir condition. These 
calculated curves do not contain an estimate of the 
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effects of the wall boundary layer. The vertical bars 
through the data points represent the variation with 
time observed during the individual experiments. At 
the lower temperature of 1,400°K, where the real gas 
effects are negligible, the observed static pressure ratios 
along the axis agreed reasonably well with the perfect 
gas curve. Fora reservoir temperature 70 of 4,000°K 
and a pressure of approximately 500 psia, the static 
pressure ratios agreed well with the values calculated 
for the case of the air in equilibrium throughout the 
expansion process. Thus, for this particular initial 
condition, the recombination was sufficiently rapid that 
nearly equilibrium flow was maintained in the nozzle. 
The experimental results of the static pressure ratios 
measured at a nozzle area ratio of 144 as a function 
of the reservoir temperature for different reservoir 
pressures are shown in Figs. 2a and 2b. The calculated 
equilibrium and frozen-expansion results are also shown 
for each case. It is seen that in all cases investigated, 
the expansion is in equilibrium for a pressure of 1,000 
psia. Fora pressure of 500 psia, the expansion process 
is very close to equilibrium up to a reservoir tempera- 
ture of about 4,500°K, before the flow process deviates 
toward the frozen condition as the temperature is in- 
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creased further (Fig. 2a). At a pressure of 200 psia, 
the flow is roughly halfway between the calculated 
equilibrium and frozen values, as indicated in Fig. 2b. 
For this particular condition the recombination rates 
for the air constituents during the expansion from the 
reservoir to an area ratio of 144 were not great enough 
to maintain equilibrium conditions in the nozzle 
Hence some of the air constituents were frozen, so that 
the static pressure is less than that for the equilibrium 
case. 

An even greater departure from the equilibrium flow 
was observed for the lowest pressure of 100 psia, as 
shown in Fig. 2b. At the higher temperatures for this 
pressure the static pressure ratios approached the 
frozen-expansion values. The degree of real gas effects 
i.e., dissociation, chemical reaction, and ionization—is 
the greatest for this low reservoir pressure and high 
temperatures. As the air plasma is expanded in the 
nozzle, most of the air constituents present in the 
reservoir at the entrance to the nozzle are ‘‘frozen’’ be- 
cause of the low densities. The collision frequencies or 
the recombination rates during the expansion process 
are not high enough to maintain the gas in equilibrium 
as the particles are accelerated and the pressure and 
temperature are decreased. The large difference in the 
static pressure ratios between the equilibrium and fro- 
zen expansion processes at low reservoir pressures and 
high temperatures is apparent from Fig. 2b. 

To investigate the expansion process closer to the 
nozzle throat, the static pressure was measured at an 
area ratio of 56.2 for a reservoir pressure of 100 psia 
and various reservoir temperatures. The static pres- 
sure results for this location are presented in Fig. 3. 
Here it can be seen that the flow expansion is close to 
equilibrium up to a temperature of about 3,000°K. 
For higher temperatures the static pressure results 
indicate that the expansion process again approaches the 
At this particular area ratio, the 


frozen flow values. 
departure from the equilibrium process occurred at a 
higher temperature than for the larger area ratio with 
the same reservoir pressure, as indicated by Figs. 2b and 
3. For example, at a temperature of 3,000°K and 100 
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HYPERSONIC NOZZLE 


psia the data for an area ratio of 144 lies at approxi- 
mately 44 percent of the distance between the frozen 
and equilibrium curves, while at the 56.2 area ratio it 
isnearly 71 percent. 

An axial survey of the static pressure was made for a 
reservoir temperature of 3,000°K and pressure of 100 
psia to determine the characteristic of the transition 
from equilibrium to partly frozen flow in the nozzle. 
The static pressure ratio is normalized with respect to 
the equilibrium case, as shown in Fig. 4. By this method 
the departure from the equilibrium expansion becomes 
more apparent. Once the deviation from equilibrium 
begins, the gas tends to “‘freeze’’ as originally predicted 
in reference 4. Further investigations must be con- 
ducted to determine whether a gas can be in equilib- 
rium at the test section after being partly frozen in 
the nozzle. 

For a two-dimensional wedge, the shock-wave angles 
are calculated to be quite different from the two ex- 
tremes of equilibrium and completely frozen nozzle 
expansion. Accordingly, a number of experiments 
were performed in which measurements of the shock- 
wave angle were made for several wedge angles. For 
high reservoir pressures the wedge shock-wave angles 
had values which indicated that the expansion of the 
air in the nozzle was essentially in equilibrium, as dis- 


cussed in reference 12. 


Conclusions 


The static pressure measurements in the hypersonic 
nozzle are very sensitive for the determination of the 
state of the gas during the expansion process. Static 
pressure probes with piezoelectric gages were developed 
to measure the pressure on the axis of the nozzle. 

For reservoir pressures greater than 500 psia, the 
expansion of the air in the nozzle is essentially in equi- 
librium up to reservoir temperatures of about 4,500°K. 
For temperatures greater than 4,000°K and pressures 
lower than 100 psia, the"expanding air is almost frozen. 
Thus, the operating regimes for equilibrium flow for a 
hypersonic shock tunnel used as an aerodynamic test 
device can be determined by this method. 

At a given area ratio for the nozzle and reservoir 
pressure, the process remains in equilibrium up to a 
certain temperature and beyond this temperature the 
flow process approaches the frozen case. The devia- 
tion from the equilibrium process is rather rapid with 
increasing temperature. 

An attempt is being made to relate the static pressure 
results to the detailed chemical kinetics of atom recom- 
bination, chemical reactions, and vibrational relaxa- 
tion with the use of an IBM 704 computer. 

At present a program is underway to make measure- 
ments of the same type in pure nitrogen; also, inves- 
tigations will be made to determine the effects of water 
vapor on the expansion process of dissociated nitrogen. 
It is hoped that eventually sufficient knowledge will be 
available for predicting the state of the gas in any 
given expansion nozzle or in the flow expansion process 
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10 100 1000 
GEOMETRIC AREA RATIO, A/A* 
Fic. 4. Variation of normalized static pressure ratio with area 


ratio for a reservoir temperature of 3,000°K at 100 psia 


over the surfaces of the long-range ballistic vehicles. 
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Space 


Svstem Specific Impulse’ 


FREDERICK W. ROSS* 


General Atomic Division of General Dynamics Corporation 


Summary 


The high specific impulse (/,,,) potentialities of nuclear propul- 
sive devices, which appear so attractive for space-vehicle applica- 
are known to be partly offset by the relatively 
A single measure for rating any propulsive 


tion, massive 
equipment required. 
system, based on its dynamic effectiveness and including this 
mass effect, is derived here as a system specific impulse, /,, 

I, is determined as the equivalent /,, of an idealized system 
having the same net effectiveness as the actual system but with 
equipment payload, 
(the 


zero m, (mass of all structural except 


From the expression of /,, in terms of actual /,,, 7 


my) 
ratio m,/mo, where mp is the initial gross mass), and the equiva- 
lent or ideal velocity v,, a number of important results follow: 

0.05, the ratio /,./Igp is 


ranging 


(1) For chemical systems with + 
0.9 for a satellite mission. For nuclear systems with 
up to 0.5, Ix./Isp reduces to as low as 0.3 

(2) For increasingly higher /,,, with +, and 7, constant, /,. 
increases rapidly at first, then tends toward a limit —v,/[g log 
(1 — rs)| for Tap — which shows that improvements in /,, 
are highly fruitful up to a point and that littie is gained by fur 
ther increase of /,, unless 7, can be reduced. 

(3) With representation of v, and log G as coordinates where 
G mo/my, the design point for any system has vector properties 
with /,, proportional to the vector argument. The /,, of a multi 
stage vehicle is proportional to the argument of the vector sum 
of the separate stages 

(4) For proportionate multistaging, /,, for 
identical with that for each stage, showing the vehicle propulsive 


the vehicle is 
effectiveness is not improved by staging. 

(5) The magnitude of G, (per stage) is determined mainly by 
the logarithmic form of the basic equation rather than by /,,, so 
the favorable payload-carrying capability of any stage, whether 
chemically or nuclearly propelled, is approximately the same 
(G, = 3 to 6) 

(6) Finally, the principal improvement of the higher /,,, of 
nuclear systems lies in the considerably higher maximum v, per 
stage by which more extensive missions (v7, up to 120,000 ft 
sec) can be accomplished with one or two stages rather than 
eight or ten, thus reducing the vehicle mass ratio G from G8 or 
G," to Gs or G2, where G, is approximately constant for all types 


Symbols 
| space-system specific impulse 
| a engine-propellant specific impulse 
| ss = equivalent useful impulse 
I, = applied impulse 
my, = total launch mass 
m = mass of payload 
wi, = mass of propellant 
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m My — Mp, — mh 
m, mass rate of propellant ejection 
t time 
t time to cutoff (of m,) 
Va actual velocity 
reference or ideal velocity 
Av, reference velocity per stage 
'e “effective” exhaust velocity of propellant ejection 
g acceleration of gravity 
a; resultant acceleration induced by gravity, drag, etc 
expressed as function of time 
my, / mM propellant mass ratio 
ry m/my = payload mass ratio 
r m,/m, = 1— fr — structural mass ratio 
6 m,/(ms + mM, r./(l — structural factor 
G mio / my l/ry 
r i', /Ve t,/2] sp mission constant 
n number of stages 
‘i thrust 
Subscripts 
0 overall vehicle 
/ ideal system 
s vehicle stage 
Introduction 


HE SPECIFIC IMPULSE, /,,, potentialities of nuclear 
"lh canals devices, which range from S00 sec to 
15,000 sec and over,'* make such devices highly at- 
tractive for space application. It is well known, how- 
ever, that they are massive and consequently not as ef- 
fective as indicated by their attractive /,, values. Yet 
in many analyses, as Sutton! or Ehricke,” systems are 
rated on the basis of /,, without reference to these 
mass differences. 

/,», being the total impulse per unit mass of ejected 
propellant, is an important measure of the capability of 
an engine-propellant combination to produce ‘“‘static”’ 
thrust. In any application, however, the problem is 
dynamic, the essential goal being to accelerate a pay- 
load to a desired velocity, the thrust being only one of 
A single 

I.» has 


been used as a convenient approximation for compari- 


several factors determining this acceleration. * 
measure for this dynamic condition is needed. 


sons but is so usable only as long as the masses of pro- 
pulsive equipment and structure such as rocket engine, 
pumps, or tanks (needed to utilize the propellant) are 
sufficiently small and similar—for example, in a com- 
two systems utilizing nearly equivalent 
In evaluations of types 


parison of 
liquid chemical propellants. 
with widely differing relative masses, however, /,, alone 
can be quite misleading. Examples include: 

(1) Comparison of a liquid hydrogen with a kero- 
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Fic. 1 Variation of I,./1 


sene system. The liquid hydrogen has about one 
twelfth the density of kerosene and hence requires extra 
mass for the much larger tank volume as well as 
special equipment to handle the low temperatures and 
evaporation. In any application, these added masses 
are necessary and must also be accelerated before the 
36 percent higher /,, can be utilized. 

Evaluation of the comparative merits of a nu- 
clear propulsive system. As the 


equipment needed for the nuclear system is relatively 


(2) 
mentioned above, 
so much more massive than that for the chemical sys- 
tem that the propellant /,, alone ioses all meaning as a 


measure of the system capability. 


(1) Single Stage Vehicles 


1.1) A Measure for Comparison of Propulsive Systems 


In this paper, a single measure is derived for compar- 
ing such widely differing schemes of propulsion by de- 
fining a vehicle or system specific impulse (repre- 
sented by /,,) based on the ability of the system to ac- 
complish a reference mission. /,, has the advantage of 
being a single measure and is derived in a form to be 
used as an equivalent /,,. 

It is observed first that /,, is a measure of the system 
input or applied impulse, and that /,, is a corresponding 
measure of the input of an ideal system which has the 
maximum capability for accomplishing the reference 


mission. 


1.2) Useful and Applied Impulse 

The end purpose of any reaction propulsive system 
is to impart an equivalent useful impulse, /, r, to 
is the ideal or 


= my 
a useful mass or payload, m), where 2, 
reference velocity defined later. 

To obtain /,, two conditions must be met: 


(1) an 


iunpulse of magnitude 


for typical values of 


839 








MASS 











structural mass factor and stage mass ratio 


I, = Sf m,v, dt 


must be applied to the system by an 
a portion of /, must be used in accelerating the ‘“‘non- 


MW, 


“engine’’; (2) 


useful’ mass of “‘engine,’’ propellant, and structure 
necessary to operate the engine and “‘carry”’ the useful 
load. 

The specific impulse being defined® as /,, v,. £, we 


then have J, = (gl,p)mp. 

/,», then, is a direct measure of the capability of the 
propulsive system to produce /, rather than /,. Our 
purpose here is to determine a measure for the capability 
to produce /,. This is the capability of the propulsive 
system (including necessary structural mass, etc.) to 


accomplish the mission. 


1.3) Reference or Ideal Velocity 


The equation of motion of a_reaction-propelled 
vehicle utilizing propellant mass castoff at a rate of 


m, may be written 


(Mio Mt yt) [(dv,/dt) + a;] Ni pd, (1) 
We define the reference or ideal velocity as 
+7 
v; = Va + a; dt 
/ 0 
27 - 
Mpd, dt m 
= . v, log 2) 
JO My — Ml m, +- mM 


Beyond its role in the definition of v, as in Eq. (2), % ap- 


pears no further in the subsequent theory. The velocity 


1 
loss term, f a; dt, has been evaluated by A. F. Don- 
0 


ovan,* for example, for typical high-thrust systems 
(7'/m) > g); and 
evaluations 


ground launch an 
analysis and summary of 
low-thrust systems operating from near-orbital tra- 


jectories with 7m as low as 0.0001g, are given by 


capable of 


typical for 
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Fic. 2a Reduction of /,, introduced by r, for v; = 12,000 ft/sec. 


J. H. Irving.‘ Interestingly, the term tends toward 
zero as 1'/mo approaches zero. 


Eq. (2) gives 
I, = Iq(m/my) log mo/(mm + ms) (3) 
which brings out the distinction between J, and /,. 
Since v, = /,/m,, the reference velocity is the useful 
impulse per unit mass of useful load and we concentrate 
our analysis in terms of v,, which from Eq. (2) can be 


written 
Vv. = —U, log (Y; + rs) (4) 
= —v, log [6 + r,(1 — 6)] (5) 


1.4) Equivalent Ideal Propulsive System 


From the fundamental nature of a reaction propulsive 
system, there would be no thrust unless a mass of 
propellant m, is carried aboard and cast off at a rate 
m, The most ideal condition possible then is when 
only m, and m, are present and m, and hence r, or 6 is 


zero. Eqs. (4) and (5) then become 
v, = —v,; log r, (6) 
= v,; log G (7) 
where G = 1/r,. 
ms; 
AGE V,= 36 
* 3 
t =( 4 
4 
. 1,70-5 _ | 
(f ; Ge 
fi 2 a 
F asia | ~~ 2000 3000 4000 5000 6000 
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Fic. 2b. Reduction of /,, introduced by r, for v, = 36,000 ft/sec. 
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(1.5) System Specific Impulse I,, 

From Eqs. (5) and (6) we have 

ve log [6 + 7(1 — 4) 
log 7; 

Thus, v,; determines the magnitude of an equivalent 
propellant exhaust velocity which will impart the same 
useful J, (= my,) to the equivalent ideal system 
(6 = () as is imparted to the actual system by 2. 

Since /,, = v./g, we may define the system specific 
impulse as 


i Vei/ £ (9) 


Eqs. (7) and (8) give, respectively, 


I,s = v,/(g log G) (10) 
I. = I,»} log {6 + ri(1 — 6)] /log ni} (11) 


It is easily seen that J,,/J., > 0 as r,;—> 0. A graph of 
I,;/1sp) as a function of G is given in Fig. | for repre- 


sentative values of 6. 


OC 
re m/m, 

3000 

ss 

2000 7 

000 Cy 

3 ~ 
Fic. 2c. Reduction of /,. introduced by 7, for », 120,000 
ft/sec 
If we define a ‘“‘mission constant’’ as c = 2,/, 


?,/gI sp, then Eq. (4) gives 


rte = 1 — tr = 6 + ri(l — 3) e” tty 


and Eq. (10) becomes 


I,, = —v,/[g log (e ‘ — r,)] (13) 


(1.6) Limiting Characteristics of I,, 
We examine two special cases of Eq. (13): 


(1) Isp > ©. Since e * —> 1, /,, asymptotically 


approaches a maximum 


(Is)mar = —2,/[g log (1 — r,)] (14) 


= —yv,/g log [1 — 6(1 — 7,)] (15) 


This shows that /,, is limited to a maximum by 7,, 
the structural mass ratio (or by 6), no matter how high 
hegit 

(2) When the term (e ‘ — r,) of Eq. (13) approaches 
zero, /,, > 0 and we have a minimum cutoff 
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(Isp)e = —0,/(g log rs) (16) 


I,, has practical value only when /,,> (Jsp)-. 
The several curves given in Figs. 2a, 2b, and 2c are 
based on Eq. (13) and illustrate these two limiting con- 


ditions. 


(1.7) Physical and Geometric Interpretation of I.,, 


In addition to its physical meaning of net propulsive 
effectiveness of the system, /,, has a fundamental geo- 
metrical interpretation. 

Since the system net output—i.e., J, = mw,—is de- 
termined by v, [see Eq. (4)] and since the ratio of gross 
mass to payload is one of the most important design 
parameters, the functional relation of v, = gl,, log G 
defined by the fundamental Eq. (10) is basic. 

Represented in a two-coordinate system with 
and log G as coordinates or vector components, any 
design point may then be considered as a vector. 
gl,, is the slope of a line through the design point and 
the origin and hence is the argument of the vector. 
This fundamental characteristic of J,,, which is illus- 
trated in Figs. 3 and 4, will be treated more fully in 


Uv, 


Section 2. 


(2) Multistage Vehicles 
(2.1) Limitations of Single Stage Vehicles 


For constant J,,, from Eq. (5) it is seen that 2, 


increases as r; decreases reaching an asymptotic limit, 
for 7; = 0, of 

(U,)maz = £1 sp log (1/6) (17) 
If the desired v, is greater than this, multistaging is 
necessary. 

Furthermore, from Eq. (11) we observe that as 7, > 0 
and v, > (%;) maz, Iss > 0, as indicated in Fig. 1. Conse- 
quently, if the system is to have a reasonable magnitude 
of effective capability—i.e., of J,,, it must be operated 
well below (%,)maz, thus requiring a number of stages 
somewhat greater than v,/(¥;) maz. 


(2.2) Analysis of I,, for Staged Vehicles 

With Av, as the increment in v, per stage, for n 
proportionately equal stages of the Malina-Summerfield 
type,® Eq. (5) gives 


v, = ndv, = —ngl,, log [6 + ms(1 — 6)] (18) 
= gl,,(n log G,) (19) 
where G, = 1/r., in which r7;, is now the ratio of stage 


payload, m,, the mass of all stages it “‘carries,’’ to 
Mos, its total mass at launch or stage separation in- 
cluding ms. 

Since for proportionate staging the ratio of actual 
payload to launch mass is 


ro = (ris)” = 1/G (20) 
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Eq. (19) becomes 
v, = gl,, log G (21) 
and Eq. (11) gives 
nI ep log [6 + ris(1 — 4)] 
log ("s)” 
log: [6 + Vis\ l - 6) 


«i (22) 
log Tis 


Is = 


Eq. (22) is identical with Eq. (11) and shows /,, to be 
independent of the number of stages. 
Using Eq. (12) we also have from Eqs. (21) and (22) 
v, 
i re (23 
tg log [(1 — 4)/(e im — §)]"} 4 
Ju(c/n) 


Ie = ———— 
{log [(1 — 6)/(e~“" — 6)]} 


[- 1 = 300 
| 8 = 0.07 


| 4 
[LOG G, = 4106 G, =106(6,/"—— xe 
| | 











LOG G 
G=GROSS MASS 
PAYLOAD MASS 


Fic. 4. Vector summation for proportionate and nonpropor- 


tionate staging. 


where now c = nAso that Eq. (24) is also independent 
of m in agreement with Eq. (22). 


(2.3) Geometric Interpretation of I,, for Staging 

From Eq. (19) and its graph given in Fig. 4, we ob- 
serve that staging constitutes a ‘‘stretching’’ of both 
abscissa and ordinate by the factor . Although higher 
magnitudes of v, are attainable by staging, this is 
accomplished by increasing the magnitude of G and 
hence by reducing or downgrading the payload ratio 
(= 1/G), by the exponential power of n (to 1/r;,”). 

At the same time the slope of the line through all 
points A, B, C, and D (Fig. 4) is still the same. Since 
this equals g/,,, this illustrates the conclusion that 
staging does not affect the propulsive effectiveness 
Tas. 

If in Eq. (21) we substitute 1/r% = (1/r,,)” for G, we 
have 


v, = gl. log 1/rm = ngl.s log 1/1), 
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TABLE 1 
Minimum Number of Stages 














n 
c= 

6 log 1/6 0.5 1.0 2.0 5.0 10.0 
0.05 3.00 1 1 1 2 4 
0.1 2.30 1 1 l 3 D 
0.3 1.61 1 1 2 4 7 
0.3 21 1 1 2 5 9 
0.4 0.917 l 2 3 6 11 
0.5 0.694 1 2 4 8 15 
0.6 0.511 1 2 4 10 20 





apparently showing the z-stage vehicle to have an 
I,, n times that of the individual stage. This appears 
to contradict the conclusion obtained from Eq. (22). 
However, v; = nAv, as noted previously, so the two 
conclusions are consistent and /,, for an m-staged 
vehicle is independent of the number of stages. 

If the staging is nonproportionate, instead of an 
n-factor stretching we have 


n 
v7, = , Avs; 
sad 


and the resultant of a typical design point such as D on 
Fig. 4 can be represented as a resultant vector with 
(g/ss) p as the argument or slope given by 


n 


In = 2, A0n/8 log G, 
s=0 


0 


with 
log G = > log G, 
s=0 


(2.4) Indirect Influence of Staging on I,, 

If 6, r;, and J,, are allowed to change with 2, it is 
seen from Eq. (23) and the curves for = 1, 2, 3, ete., 
in Fig. 4 that designs can be selected to have /;, greater 
(see point A’). 

In the limit, with Av, > 0 andu—> o~, 
increased to a maximum, the tangent at the origin, for 
which it can be shown that, from Eq. (24), 


gl,, can be 


I,5 > (1 — 5)Isp (25) 


(2.5) Minimum Number of Stages 
For I,, to be positive, the term e “” — 6 in Eq. (24) 


must be positive so that only cases having 





TABLE 2 
Typical Magnitudes of c (= 7;/gIsp) 





c for Various v; (ft/sec) 








Isp 30,000 60,000 90,000 130,000 200,000 
400 2.33 4.65 7.00 10.1 15.5 

600 1.55 3.10 4.65 6.70 10.3 

800 : es 2.33 3.50 5.05 7.70 
1000 0.93 1.86 2.80 4.03 6.22 
2000 0.47 0.93 1.40 2.02 3.11 
4000 0.23 0.47 0.70 1.01 1.55 
6000 0.16 0.31 0.47 0.67 1.03 
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are of practicalinterest. This determines the minimum 
number of stages required for any mission constant c. 
Table 1 lists these, and Table 2 lists magnitudes of c in 
terms of v, and J,». 


2.6) Preferred Number of Stages 

As c gets larger, the minimum number of stages as 
determined by Eq. (26) increases, and for any number 
n above that minimum the ratio J,,/J;, approaches 
rapidly to an asymptotic value. Since it is in the 
interest of low cost, less complication, and higher 
reliability to keep m low, then m must be confined be- 
tween these rather narrow limits. Three typical ex- 
amples are illustrated in Fig. 5. 


(2.7) Minimum Gross Mass Ratio 
From Eqs. (21) and (25) the limiting gross mass, 
when 7—> © is 


, Vr/olss ssai—é 97 
Go. =e —— tall (27) 
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For the preferred number of stages 
. Wa. - c/(1—8) ; 
Ge = teGe == Jee ( 5 (28) 
where 


[1 — 6)/(e~”" — 6)]” 

fe = “= (29) 
e 

The dependence of /, on 7 is illustrated in Fig. 6 for two 


long-range space missions. 


Conclusion 


Introduction of the concept of J,; as a measure of the 
dynamic effectiveness of a space system makes it 
possible to rate systems of various types of propulsion 
with different J,,, 6, v;, and 7, on a single common scale 
for direct comparison. 

Furthermore, (1) it separates the system net ef- 
fectiveness J,, of a multistage vehicle from the down- 
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Mass-increase factor for n < 


grading of payload-carrying capabilities introduced by 
the staging—.e., from the exponential reduction of rj; 
to r,;"—and (2) it shows that the system effectiveness 
for each stage is equal to that for the overall system if 
staging is proportionate, or to a ratio of simple sum- 
mations if the staging is nonproportionate. 

Thus space vehicle analysis involving staging can be 
more directly concentrated on the characteristics of 
the individual stage with little more than adding or 
“stacking” of the stages to obtain the necessary 
v, = nAdv, where Av, is the limited reference velocity in- 
crement per stage. 

Following this procedure, it is observed from Fig. 3 
that for all curves for single stage vehicles—i.e., with 
n = 1, the “knee,” of the logarithmic curves all occur 
at approximately the same value of G This 
effect is brought out clearly by representing v, non- 
dimensionally by combining Eqs. (5) and (17) as 


1.e., G.. 


Vr log [6 + G,-'(1 — 6)] 


log 6 


(v,) mar 


This equation is independent of /,, and v, and shows 
that the relative velocity increment per stage is a func- 
tion only of 6 and G,. A graph of Eq. (30) given in 
Fig. 7 further shows the rather limited effect introduced 
by 6 when 2,/(%;)maz is considered as a function of G,, 


(Continued on page 854) 
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Hypersonic Flow Over an Elliptic Cone: 
Theory and Experiment! 


ROBERT L. CHAPKIS* 


Guggenheim Aeronautical Laboratory, California Institute of Technology 


Summary 


By applying hypersonic approximations to Ferri’s linearized 
characteristics method, simple results were obtained for the 
shock shape and surface pressure distribution for an unyawed 
conical body of arbitrary cross section. Calculations were 
carried out for an elliptic cone having a ratio of major to minor 
axes of 2:1, and a semivertex angle of about 12° in the meridian 
plane containing the major axis. An experimental investigation 
of the flow over this body conducted at a Mach number of 5.8 
in the GALCIT hypersonic wind tunnel showed that the surface 
pressure distribution at zero angle of attack agreed quite closely 
with the theoretical prediction. On the other hand, the simple 
Newtonian approximation predicts pressures that are too low. 

Surface pressure distributions and schlieren photographs of 
the shock shape were obtained at angles of attack up to 14° at 
zero yaw, and at angles of yaw up to 10° at zero pitch. At the 
higher angles of attack the Newtonian approximation for the 


surface pressures is quite accurate. 


Symbols 
a = speed of sound, referred to limiting velocity 
Cp = pressure coefficient 
M, = free-stream Mach number 
R = gas constant 
Ss = entropy 
u,v,w = velocity components in spherical coordinates, re- 
ferred to limiting velocity (see Fig. 1) 
V = local resultant velocity, referred to limiting velocity 
Vi = undisturbed velocity, referred to limiting velocity 
r,0,@ = spherical coordinates (see Fig. 1) 
a = angle of attack 
¥ = ratio of specific heats 
¥ = angle of yaw 
Subscripts 
n,m = indices of summation 
0 = properties of the basic axisymmetric conical flow 
c = properties at the body surface 
Ss = properties at the shock wave 
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(I) Introduction 


S TUDIES OF THE FLOW past elliptic cones are of inter- 
est for two reasons: (1) as Van Dyke has pointed out,! 
the elliptic cone will probably become a standard of 
comparison for supersonic flow past bodies without 
axial symmetry, just as the circular cone is used as a 
of comparison for supersonic flow past 
(2) elliptic cones may have 


standard 
bodies of revolution; 
important aerodynamic advantages over circular cones. 
For example, both theoretical and experimental 
investigations? * have shown that an elliptic cone may 
have significantly higher lift/drag ratios than a circular 
cone of the same cross-sectional area per unit length. 

Most of the theoretical investigations of the flow 
about elliptic cones are based upon linearizing as- 
sumptions which are not valid at very high Mach 
numbers. Even Van Dyke’s second-order theory,! 
which proceeds from slender body theory and includes 
the effect of the leading nonlinear terms in the equations 
of motion, cannot be expected to give good results in 
the hypersonic speed range. f 

At hypersonic speeds the well-known Newtonian 
approximation (or a suitable modification®) has been 
quite successful in predicting surface pressure dis- 
tributions, provided that the component of Mach 
number normal to the surface is of order unity, or 
larger. But for ‘‘flat’’ bodies, such as a delta wing 
or an elliptic cone of significant eccentricity, this 
restriction means that the free-stream Mach number 
must be extremely high, or the angle of pitch (or yaw) 
must be large. For example, the Newtonian approxi- 
mation states that the shock wave coincides with the 
body surface, at least to the first order. However, on 
physical grounds the cross section of the shock surface 
for an unyawed elliptic cone is expected to have a 
smaller eccentricity than the body cross section, 
except possibly in the limiting case (y — 1)M/ 
This fact must have a significant effect on the surface 


2—> ~, 


pressure distribution. 

In view of these criticisms, it seems desirable to work 
out a solution for the unyawed elliptic cone at hyper- 
sonic speeds directly from the gasdynamic equations 
of motion. One attractive approach to this problem 
is the ‘linearized characteristics method’’ developed 


- t It turns out that Van Dyke’s theory does give good results, 
at least for the elliptic cone which was tested at M = 5.8 (see 
Fig. 16). 
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by Ferri and co-workers,” * which considers the flow Vv 
about a body to be a perturbation of a known basic 
nonlinear flow field. Thus, the flow about an elliptic U 
cone is considered as being a perturbation of the 
known flow field about a circular cone. Only linear 
terms in the perturbation quantities are retained in t a 
the differential equations and boundary conditions. 
It should be noted that the linearization is with respect Pp 
to deviations from a known basic flow field which is 2 
“close’”’ to the actual flow field, and not with respect 0 
to deviations from the uniform flow upstream of the 
shock wave. This method should be applicable for 
hypersonic as well as supersonic flows, provided that 
the ‘‘exact’’ basic flow field is known. 

Ness and Kaplita’ utilized Ferri’s scheme, and the 
results of Kopal’s computations for the unyawed Fic. 1. Coordinate system and velocity components. 
circular cone, to calculate the perturbation velocities 











yf inter- 
ed out,! 
dard of 


rp. for any conical body which employs the circular cone simplified. In fact simple algebraic expressions are 
; flow field as a basic flow. Their calculation entails obtained for these perturbation velocities (Section II). 
; ee a step-by-step numerical integration of the governing Once these quantities are known the surface pressure 
scales differential equations for each particular case. Now distribution and the Fourier coefficients for the shock 
sesiaiital for hypersonic flow this difficulty can be circumvented shape are readily calculated. 
arg: by employing the hypersonic approximation for the The ultimate test of the theoretical analysis must 
rae: flow over a circular cone obtained by Lees.* By be the comparison of it with an “‘exact’’ solution or 
oats expanding the velocity components in a Taylor's with experimental results. Several experimental in- 
ie flow | tes in the conical ray angle, Lees obtained simple vestigations have been made of the flow over elliptic 
we al approximate expressions for the velocity components, cones at low and moderate supersonic Mach numbers 
etal shock wave angle, and pressure coefficient. These (up to about Mach 3). %, 10 However, no experi- 
heated approximate expressions give results which agree very mental results were available at hypersonic speeds. * 
ied well with the actual values computed by Kopal, Therefore, an experimental investigation of the flow 
Nee provided that the hypersonic similarity parameter over an elliptic cone was carried out in the GALCIT 
omar: K = M6, is greater than about one. By substituting M = 5.8 wind tunnel, in order to obtain surface pres- 
sate. these approximate expressions for the circular-cone sure distributions and shock wave shapes at zero angle 
rer velocity components into the differential equation of attack and also at various angles of attack and of 
mee governing the perturbation velocities for a body of yaw. The description and results of the experimental 
re yo noncircular cross section, this equation is greatly investigation are presented in Section III. 
Mach 
ity, or 
per . (II) Theoretical Investigation 
, 
umber | 4. Resume of the Linearized Characteristic Method 
r yaw) By combining the continuity, momentum, and energy equations for a nonviscous perfect gas, the following equa- 
proxi- | tions are obtained for the case of conical flow: 
th the : m P “i sae a a 
rer, on “ (2 - Bw) + v cot o +5 (1 -")4+>- 2 (1 -“)-" (= _ +o) —~ 0 (1) 
urface a® 06 a’® sin 6 O¢ a* a* \sin@0¢ 08 
_— a a? Os ; ow Ou Ov : . ; (2) 
ction, yR a6 = vsin 0 7] —U de —v de + wu sin 6 + wv cos 6 2 
—> OO. 
urface a* Os Ou ow w Ov ; i 
ae”6° “ES "sae "YO (3) 
. work , : — 
vy per- where the directions of u, v, w, 0, and ¢ are as shown in Fig. 1. 
ations Immediately downstream of the shock front—i.e., at @ = @,—the following relations exist: 
oblem : 
loped (u)e, = Vi cos 6, (4) 
esults, * After this paper was prepared, experimental results on two elliptic cones at M = 6 obtained at the Polytechnic Institute of 


8 (see Brooklyn Aerodynamic Laboratory were reported in references 13 and 14. These results were compared with the Newtonian ap- 


proximation, with an ‘equivalent cone’ method, and with an extension of the linearized characteristics method. 
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a. y—1/1— V,? cos? é, — V;? sin? @, sin? B aes ts 
(vy) = — ' aan aiee era ——— }] — V,; sin 0, sin? B (5) 
++ 1 V; sin 0, ’ 
seule . ¥—1i. oa si ee SE 
(w)e, = —Vi sin @, sin B cos 8B + — (1 — V,? cos? 0, — V;? sin? 6, sin? 8) ———— (6) 
7? % V; sin @, 
where 0, is the shock wave angle, and £ is defined by 
tan 8 = (1/sin 0,)(d0,/d@) (7) 
The shock wave angle @, is expanded in a Fourier series: 
0, = Bos + On. cos nd + DOOns sin mo (8) 


where 4, is the shock angle for the basic circular cone. In Ferri’s analysis only first-order terms in 0,5 and On; are 
retained in the differential equations and boundary conditions. To this approximation, series representations of 
the velocity components which are consistent with the conditions just behind the shock front are as follows: 


u=uUt+ Onsen cos nd + > Omsttm sin mo (9) 
DV = + LOnsYn COS NP + YOmsVm Sin mp (10) 
w= 18 ,5Wy sin nd + } m MO msWm cos mp 1) 


where W,, Um, ete., are functions of 6. 
By expanding uw, v, and w in Taylor’s series around @ = 4, and recognizing that 0, = 00, + (@, — Os), the shock 
conditions [Eqs. (4)-(6)] can be made to yield the following relations for u,, u», etc., at 6 = O,., independently of 


n and m: 


(thn) 3 = (Um) = — Vi Sin Bos— (V0) ms (12) 
P . ies l . l a V,? cos? Gos Ovo 
(Un) 6, == (Um) 65 = — - COS Oo5 (21 os eee — : (13) 
¥+1 Vi sin? 6; 00 / ¢,, 
(Wn) es = (Wm), = Vi + (1/sin Box) (v0) (14) 
Here mp and vp are the radial and normal velocity components corresponding to the flow about a circular cone of 
semivertex angle 0-. The shock wave angle for the flow about the circular cone is %,. At @ = Os: 
(0); = Vi cos 4, (15) 
: y¥— 1 (’ — V,? cos? *) (*) ] Uo + Uo cot | 
V0) = — eee = —| m+ _—T (16), (17 
( _ Y -4 l | 1 Sin O05 00 Pes , | (Vo ao)* Pes 7 ) ( ) 


When the series for “, v, and w are substituted into Eq. (1) and only first-order terms in @,, and 6,,, are retained, 
three independent equations result. As expected, the quantities 7) and v satisfy the differential equation governing 
axially symmetric conical flow: 


uo (2 wid (vo ao)? | Vp cot 0 + (Ov 06) [1 oe (vo/ao)?| = 0 (18) 


The second equation is 
Ov, NW» Vo \? Ov», (vo ‘ag) (uo/ao) i (vo ‘do)? cot 6 , Uo \? 
ee + Un cot 6 = ZU n og P = et . ae = 2 + 7 1) Un -+- 
06 : sin 0 ao/ 086 1 — (v/do)? \ Qo 
is is (vo/ao)(uo/ao) + (v/ao)? cot 6 
7 (Uo, ao)? = — . . 


1 — (vo /ao)? 


(y — 1)(uo/ao) (v% an) un (19) 
f 


The third equation is exactly the same as Eq. (19) but with the subscript ” replaced by m. 
In order to solve for u,, v,, and w,,, two more relations between u,, v,, and w, are needed in addition to Eq. (19). 
Ferri obtains these relations from Eqs. (2) and (3) and the momentum equation in the radial direction. To the first 


approximation in @, and @,,: 
Vn = Ou,/O8, Vo = Ouo/00 (20) 
and —(ao?/yR)s’ = v0 sin 0(0w/00) + Uolltn + Von + UW, Sin 6 + WwW, COS O (21) 
where s’ = (ds/d0,)o, is the rate of change of entropy with respect to shock angle. By utilizing Eq. (20), Eq. (21) 
can be written as 
— (ay?/yR)s’ = uo(tn + Wy sin 0) + 19(0/08) (uw, + w, sin 6) 22) 
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/ 6 y= te = 

: Ss ie. y 3 ay 

or Un + W,S5ind = — ( ) agi! 7 D(—v sin @)!/? [ - - dé (23) 
yR J 6,, Vo(—to sin @)'/? 


The relations between u,,, 7», and w, are the same as Eqs. (20) and (23) with the subscript m replaced by m. Eggs. 

(19), (20), and (23) plus a knowledge of the values of uo, vo, and 6, enable one to compute u,, “4, etc., numerically. 
The boundary conditions at the surface of the body enable one to determine the shock coefficients, Am. and Ons, 

once v, and v,, are determined. The boundary condition at the surface of the body is that the velocity component 


normal to the surface of the body be zero, or for conical bodies: 
(v/w)e. = (1/sin 6.)(d0./d¢) (24) 
The shape of the body as defined by 6, = 6,(¢) is now also expressed in a Fourier series as follows: 
0. = Doe + » » Onc COS NP + _ Ame SiN mod (25) 


One observes that (v/w)¢,, is of order 6,., and w is of order @,,, so that v is of order 6,,.0,-—1.€., (v)e, O to the first 
order. Now v and w are expanded in Taylor’s series about 0 = ®,. If only first-order terms in 4c, Amc, Ons, and 
6,,; are retained, one finds 


Ons = Znc(tUo/Vn) bo¢ Ams = 2Ome(t0/Vm) Oe (26) 


The pressure distribution is obtained from the expression 


im . (3 = ~) —As/R | 97 
>= M2 L\i — Vit et) 


The square of the magnitude of the velocity on the surface of the body is 
V? - [ (200) b¢ si Zz, Bns(Un) Oe cos np + oe Ams (Um) Oe sin mo |* + oe NI ns(Wn) bo¢ sin nop - > i MB ns(Wm) bo cos mo |? (28) 


since to the approximation accepted, (v)», = 0. Ferri shows that 0s/06 is of order (@,,)*; therefore, the entropy 
behind the shock front is taken as constant for each meridian plane but varies from plane to plane. 


(B) Hypersonic Approximation 

In this section, approximations valid for hypersonic flow are applied to Eqs. (19) and (23), which will enable simple 
algebraic expressions to be obtained for u,, v,, and wy. 

First of all, Eq. (19) will be written in a more convenient form. Solving for 0v/06 from Eq. (18), it is easily seen 


that 


Vo ( =) (vp/ao)(uo/ao) + (v/do)? cot 6 
= 1 5 ——— oi ; 
06 1— (vo/ao)? 


2 


ao 


Therefore Eq. (19) can be written in the form 


Ov, + c t 6 4 2) 4 NW (* ) E Ov, 7X Vo ( 4 mn) 2 oo ( l \ (* | a 
. UV, CO su» : = 10 4 = Un 
06 sin 0 ao/ O08 ao? 06 , * \ao 


%\? | Ov uo\/v 
v0 ev et ( 0 
+ — (w+ )ar- »( )( )| u, (29) 
| (*) Ao” O86 ao/ \do 
The coefficients of 7,, v,, and Ov,/06 on the right-hand side of Eq. (29) involve only the axisymmetric conical flow. 
This equation can be greatly simplified for the case of hypersonic flow by making use of the series representations 
for u% and v which were derived by Lees in reference 8. The series are the following: 
Uo = (to) %,[1 — (8 — Bac)? + 1/3 cot Oc — Boc)® — as(O — Boc)* +... JQ (30) 
= (10) 6. —2(8 —_ Doc) + cot 60-(0 _ Boe)? eae 4a,(0 = o-)* + 2% a f 


Oa 
_~ 


1 > Uo” 
where a, = — cot? , + ( ) 
a " " 3(y — 1) \1 — we?/ a, 


In this analysis it is assumed that (@ — 6.) cot 6. is much less than one, which implies that the shock wave is ““close”’ 
to the body surface. If this assumption is made throughout the analysis, and if Eqs. (30) are substituted into Eq. 
(29), the following equation is obtained: 


Ov, my edie & NW 
= Un CO Qu» - = 
06 sin 0 


8 Uo? Ov Uo” ; 
(0 — O.)? — (0 — Oo.)Un 2 ( ) (0 — 6 Jus | (31 
1 ( :). | te) reyi] t T 1 — 07/6, ) 


Y — l — Uo” 
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‘i ( Uo? ) nwivg-1 ao * 
owever 8 = - (52) 
? L = uo" 8c 27 +1 0. — 0 


Thus Eq. (31) reduces to the following: 


ws 
ma 


Ov, 4 ii ed nw 4 (; - *: | - 
- Vn CO Zu» — => Un CO or 
o0 sinO0@ y+ 1 \OQoc — 4 . 


In a similar manner, Eq. (23) relating u, and w, is simplified by utilizing the series for w and vp» valid for hyper- 
sonic flow. If, in addition, one makes the slender-body hypersonic approximation, then cot 6 = 1/sin@ = 1/8, 
and Eq. (23) reduces to the following form: 


. y¥—1f/1— uo? . 6 Oo. \1/? f 80 — Oo. \1/? 
Un + W,5iné@ = —- és i a ; (34) 
yi U0 Bec yR boc 6 00s = Hoe 


Now Eqs. (20) and (34) are used to write Eq. (33) in terms of w, only: 


Oru» 1 Ou, n? 4 0 — Oe \ 1 Oup, 
—— pn —— io — it, = = = es —~ 
06? 6 O06 0? y + 1 \0,, — 9c} As O08 


n? E —1 (; _ “) 5! (- yy (“)( 0 — Boe ) / (35 
— {| - - - = : 39) 
g? 2 0 %0¢ yR Hoc 6 As —_ Io 


It turns out that the second term on the right-hand side involving ds/d@, makes only a small contribution to the 
solution for “,, and the contribution of the first term involving 0u,,/00 is even smaller. Therefore, a first approxima- 
tion for u, is obtained by neglecting the first term on the right-hand side of Eq. (35). Also, it is assumed that n?/@ 
is much larger than 2, which is certainly true for a slender body. If these assumptions are made, then the solution 


of the homogeneous equation for 7, and v, is as follows: 


some = —ove[ 4) ~ (Ted mel 8) +N 
Un)homog. = — ~ (Un) Os : — ~ (Un) bs | 
‘ 2n ™ 0 os 2 ' i] os 
he [(G)"- GQ) ]-see G+ QR, 
ae Cees ae i Un) os — 7? = = ee n) bs | noe p 36 
n 2 1 \%G 8os go \ee Ho) + \a, | “7 
] ‘ Caan eye n F Mos os \” * e\r | 


The particular solution for u, can be obtained by the method of variation of parameters, as follows: 


— 1/1 — up s’ Mn : 6\” ag O0.\/2/ t — Oe \1/2 
(4,,) partic. = = ( ) < : f | (7) — (5) |: —_ ( ) (; 7 ~ ) Ja (37) 
2 Uo Boe yR 260, J 8 t 7) t As — 9b¢ 


os 


A correction to the first approximation for u, is derived by substituting the expression for 0u,,/0@ from Eq. (36) 
into the first term on the right-hand side of Eq. (35). The method of variation of parameters is then used to obtain 
a particular solution which is a correction to the first approximation for u,. This procedure is carried out in Ap- 


pendix 2. 
According to Eqs. (12), (18), (15), and (16) the values of (w#,), and (v,)@, depend only on the axisymmetric 


conical flow, as follows: 


(tn) Hs = — (Uo); tan Os — (Vo) Hs (38) 
: | el oe ov 
(Un) &s = —2 ¥ rt 1 (140) Os = (Uo) Os cot 0; — (=). (39) 
By utilizing Eqs. (17) and (31), one finds that 
(tn) 605 = — (U0) 90s (Boe Bos)? 
(40) 


( ) ~w ( ) | 4 (1 O65 a “| 4 ( ) 
ne, Gt Goda, | ——— 11 + ——— ] | & —— (ee, 
™ Ly +1 Occ yt1°°" 


When these relations are employed in Eqs. (36), the expressions for (u,), and (v,)», take the following form: 





* This relation is derived in Appendix 1. 








hyper- 
= 1/6, 


(35) 
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|. (36) 
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Bo f 2 | (Fe) (*)" , Ps n | (Fe) : 4 (*)"" ) 
n ly + l Boe Bos 2 Ay. Ao. { | 
Vn 2 &,\*T* Go. \"~? n ee Bo-\" +} 
~ + — ) 
Uo %- y¥+ 1 Bo. 60s 2 Box A. 
and from Eq. (34): 
(“) l (‘“") y—1 (? - “| (=) 
—_> — — (49 
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and K, = 1,6... This expression for 
y-—1 (’ — ~) (*,) 
2 Uy” Boe NR 
is derived in Appendix 1. 


The expressions for %,,, ¥,, and W,, are the same as those for u,, v,, and w,, but with the subscript ” replaced by m. 
By using Eq. (27) for C, and Eqs. (41), (42), and (26) the surface pressure distribution is calculated. Shock shape 


F a. 
P 
s 

Il 


(41) 


I| 


is calculated using Eqs. (8) and (26). 
Although the present calculation is carried out for a particular elliptic cone, evidently the procedure is applicable 


to a slender conical body of arbitrary cross section at hypersonic speeds. 


(III) Experimental Investigation of Hypersonic 
Flow Over an Elliptic Cone 


A) Experimental 

The experiment was conducted in the GALCIT 
5 X 5-in. hypersonic wind tunnel, which is a closed- 
return, continuously operating tunnel. All of the tests 
were made at a fixed reservoir temperature of 250°F 
and a fixed stagnation pressure of 74 psi gage, giving 
a nominal test section Mach number of 5.8 and a 
Reynolds number per inch of 2.2 X 10°. An extensive 
description of the experimental facilities is given in 





reference 11. 

The elliptic cone model was constructed of brass, 
and the length of the model from tip to base was 
43/1, in. The ratio of major to minor axes of the 


— oe 


Angle of Attack 


elliptic cross section was 2, the major axis being 1*/, in. 
78.5° long at the base, and the minor axis */, in. 
Twelve pressure orifices of 0.016 in. diameter were 
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Elliptic cone model mounting. 


Fic. 2. Orifice location and notation. Fic 3. 
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Schlieren photograph at a = 0°. 


Fic. 4. 


located around the periphery of the cone in a section 
1%/i. in. from the tip. The angular locations of the 
orifices are shown in Fig. 2. Note that there are 
four orifices which are spaced 90° apart, two along the 
ends of a major axis and two along the ends of a minor 
axis. These four orifices permitted alignment of the 
model with the flow direction. 

The surface pressures were measured by means of 
a multitube vacuum-referenced silicone oil manometer. 

As shown in Fig. 3, the model was mounted in the 
wind tunnel on a 1/2-in. diameter sting which was 
supported by two vertical struts extending through the 
top of the test section. The vertical struts were 
individually raised or lowered by external controls 
to adjust the pitch angle of the model. It can be seen 
in Fig. 3 that the rear vertical strut is not connected 
directly to the sting, there being an intermediate 
short piece of steel. This short piece of steel could 
be moved from side to side by means of an attached 
thin rod passing through a hole in the side of the wind 
tunnel. Small corrections for yaw misalignment could 
be made by pushing or pulling the wire from outside the 
wind tunnel. 

Pressure measurements were obtained for angles of 
attack of 0°, +2°, +4°, +6°, +8°, +10°, and +14°, 
and for yaw angles of +2°, +4°, +6°, +8°, and 
+10°. In order to obtain angles of attack above 





Fic. 5. Schlieren photograph at a = 8°. 
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ten degrees, it was necessary to use a sting with a ten- 
deg. bend where the sting attached to the model. 
Schlieren photographs were also taken at these angles 
of attack and yaw. The pressure measurements are 
estimated to be accurate within +1 percent, and the 
angles of attack and yaw accurate within +0.1°. 


(B) Results and Discussion 

Schlieren photographs are shown in Figs. 4-9. By 
comparing Figs. 4 and 7 for the model at zero angle of 
attack and zero yaw angle, one sees that the distance 
between the shock wave and the body is greater in the 
meridian plane containing the minor axis than in the 
meridian plane containing the major axis. Thus, the 
shock wave is not an ellipse similar to the body, but is 
“pushed in’ toward the major axis and “pulled out” 
from the minor axis. This shape is to be expected 
on physical grounds, since there is a cross-flow from 
the high pressure sides at the ends of the major axis 
to the low pressure sides at the ends of the minor axis. 
This cross-flow tends to relieve the pressure somewhat 
at the high pressure sides and raise the pressure at the 





Fic. 6. Schlieren photograph at a = 14°. 


low pressure sides, thus causing a corresponding change 
in shock wave angle. 

The surface pressure distributions are plotted in 
Figs. 10-16, and the experimental pressure distribu- 
tions are compared with the pressure distributions 
predicted by simple Newtonian theory in Figs. 13, 
14, and 15. For a = 0°, the Newtonian theory 
predicts surface pressures which are about 50 percent 
too low on the low pressure sides of the cone (Fig. 13) 
but are very close to the experimental values at the 
high pressure side. This result is expected from the 
physical argument given above—i.e., from the cross- 
flow from the high pressure to the low pressure side. 
As a increases, the Newtonian theory becomes more 
and more accurate (Figs. 13 and 14). This result is 
also expected, since at higher angles of attack the 
component of Mach number normal to the surface is 
larger, and the shock wave on the high pressure side 
is closer to the surface of the body. 

In Fig. 16 the experimenta) surface pressure dis- 
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tribution for a = 0° is compared with the predictions 
of Newtonian theory, Van Dyke’s second-order slender 
body theory, and the present theory. Both the theory 
developed here and the slender body theory give 
pressures which are close to the experimental values, 
while the Newtonian theory predicts pressures which 
are too low over most of the surface. 


(IV) Conclusions and Summary 


(1) By applying approximations valid for slender 
bodies in hypersonic flow to Ferri’s linearized 
teristics method, simple algebraic expressions 
perturbation velocities are obtained for an unyawed 
Once these 


charac- 
for the 


conical body of arbitrary cross section. 
perturbation velocities are determined, the surface 
pressure distribution and the Fourier coefficients for 
the shock shape are readily calculated. The method 
is applied to an elliptic cone having a ratio of major 





= 0°. 


Schlieren photograph at y = 


Fic 7. 


to minor axes of 2:1, and a semivertex angle of about 
12° in the meridian plane containing the major axis. 

(2) At zero yaw and pitch the experimentally meas- 
ured surface pressure distribution is predicted closely 
by the present theory and by Van Dyke’s second-order 
slender body theory. On the other hand, the simple 
Newtonian approximation predicts pressures that are 
too low. 

(3) Schlieren photographs of the 
taken at zero pitch and yaw in the planes containing 
the major and minor axes of the elliptic cross section 
show that the shock surface lies considerably farther 
away from the body near the ends of the minor axis 
than near the ends of the major axis. This behavior 
is consistent with the “relieving effect’’ predicted on 
physical grounds, and helps to explain the calculated 
and measured surface pressure distribution. 

(4) As the angle of attack is increased the measured 
surface pressure distributions agree more and more 


shock surface 


closely with the Newtonian approximation. At the 
highest angle of attack (a = 14°) the Newtonian 


approximation is quite accurate. 
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Schlieren photograph at y = 4°. 


Fic. 8. 
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Appendix 1: Hypersonic Approximations for 
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The change in entropy across a shock wave of angle 
#, is given by the following expression: 
—(a)e/R _ | (y + 1)M;? sin? 6, i ditties 
¢ —_ 7 Seu" a one ok a 
(y — 1)M;? sin? 6, + 2 


y ‘4. 1 W/(y — DD 
E sin? @, — (y — 5 | 


Therefore, for constant y and R, 
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I(y — 1)M)2 sin? , + 2] 
oe (K,? — 1)? 
~ Os [2yKs2 — (vy — LI [(y — 1)K,? + 2] 
where K, = M,6s;. From Eq. (5) of reference 8, it is 
easily seen that 
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Bo¢ 


4 Ue? _ 2M,? 


Therefore, 
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To obtain an expression for [m?/(1 — wo”) Jo, in 


terms of (6 M,-), use is made of the fact that at 


OG = Op, 
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From Eq. (39) it is seen that 
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* This relation is easily derived from Eqs. (15) and (16). 
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Therefore, 
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Appendix 2 
Correction for First Approximation for u, and v,, 


If the expression for v, = Ou,/00 from Eq. (36) is 
substituted into the first term on the right-hand side 
of Eq. (35), this term takes the following form: 
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Y + 1 \0. — Hoc Hos 09 

Ss ( Uo? ) Ou», 
- - (0 — &-) 
y¥—1\1 — u*/«, O00 


8 Uo” 0 — A, 
= - n (A@” — Be-") 
y¥— 1\I — u*/a, 6 


| 0; 
A woe 9 Hos ” (Un) Os + (Un) Os 
n 


| Bos 
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2 n 


The method of variation of parameters is now used 
to obtain a particular solution which is a correction to 
the first approximation for u, and v,. The corrected 
expressions for u, and v, are given by the expressions 
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The values of C and D obtained by utilizing the 
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boundary conditions on u,, v,, and w, at 9 = Os, are 


as follows: 
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By evaluating the integrals in these expressions the following relations are finally obtained for C and D: 
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(Continued from page 843) 





even for a rather widespread range of values of 6. 

Thus it becomes evident that regardless of magnitude 
of J;,, whether as low as 250 sec, as for some chemical 
systems, or as high as 15,000 for a nuclear system. 
because of the logarithmic nature of the basic dynamic 
equation of an accelerating propellant reaction system, 
the payload-carrying capability (7, = G,~') is limited, 
for reasonably practical magnitudes of /,, (as indicated 
Fig. 1), to a narrow range of values approximately 
given by the inequality 


3<G,<6 (31) 


Recalling that for a vehicle of m stages the overall 


mass-to-payload ratio, by Eq. (20), is G = (G,)”, 
the inequality (31) becomes 
Fr <G¢<¢ (32) 


From this it is evident wherein the main advantage 
of massive, high J;, systems, such as potential nuclear 
types, lies; and we reach the important conclusion 
that while nuclear systems will have about the same 
magnitude for the stage mass factor G,;, with them it is 


possible to reduce the number of stages by an order of 
magnitude and hence the overall mass ratio G by an 
exponential order of magnitude. Expressed numer- 
ically for a more extensive mission (v, = 120,000 
it/sec), this typically reduces G (if G, = 4) from 4° 
to 4. This is a reduction in the gross take-off mass, 
for a particular payload mass, of more than 10°. 
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The Curtain Jet 


FREDRIC F. EHRICH* 
General Electric Company 


Summary 


A detailed analytic study is made of the curtain jet, the two- 
dimensional fluid wall used to contain support pressure on the 
underside of ground effect machines. Two variations of the jet 
are studied in detail—the bifurcated jet, in which a portion of the 
flow streams into the support pressure region, and the deflected 
jet, in which none of the flow penetrates into the support pressure 
region. Kirchhoff-Helmholtz free streamline analysis is used 
to construct the flow field, and quantitative results are presented 
for the effect of nozzle inclination and detailed geometry on flow 


requirements and support pressure differential at varying 
altitudes. 
Symbols 
a ye 2 } = identification of significant mapping points 
Co28 Os (Fig. 2) 
E, 2’, fF, } 
a = location of the nozzle source in the ¢ plane 
C = Schwarz-Christoffel mapping function con- 
stant 
G. 8 £ = integrals 
m = substitute variable of integration 
Po = stagnation pressure in nozzle 
py ambient pressure 
pe = support pressure 
Q = log of the complex conjugate velocity 
q = fluid velocity 
1, Ti, 2 = location of nozzle edges in ¢ plane 
s = location of the stagnation point in the ¢ 
plane 
t = width of the nozzle 
th, te = width of the efflux streams 
U;, Ur = fluid velocity of the efflux streams 
U = fluid velocity in the nozzle 
u = fluid horizontal velocity component 
v = fluid vertical velocity component 
W = nozzle flow rate 
Wmaz = nozzle flow rate away from ground effect 
x = horizontal physical coordinate 
y = vertical physical coordinate 
Vi, Ve = altitude of the nozzle edges 
z = complex physical coordinate 
a = inclination of the nozzle 
6 = dimensionless index of relative nozzle edge 
location, Eq. (1) 
c = complex half-plane mapping function 
n = imaginary part, complex half-plane map- 
ping function 
9 = inclination of a streamline 
y = complex conjugate velocity 
é real part, complex half-plane mapping 
function 
cy = potential function 
y = stream function 
w = complex potential function 
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Subscripts 
1 = referring to conditions on the ambient side 
of the curtain jet 
2 = referring to conditions on the support pres- 


sure side of the curtain jet 


Introduction 


| Big ALTITUDE airborne vehicles, termed ground 
effect machines, depend for support on develop- 
ment of a static pressure greater than ambient on the 
underside of the craft. In one type of design, Fig. 1, 
the pressure is contained by a curtain jet, a wall of 
fluid around the periphery of the craft, which im- 
pinges on the ground. Analyses of this curtain, or 
annular jet, have been given by Chaplin,’ Pinnes,? 
Chaplin and Stephenson,* and Strand.‘ 

Chaplin presents a very simple approximation based 
on the assumption that the jet is very thin, and of 
circular are cross section. Some refinement is added 
by Pinnes who considers a finite jet thickness, hypothe- 
sizing a free vortex velocity distribution across its 
width. Chaplin and Stephenson further elaborate 
the analysis in their “‘thick-jet’’ treatment, by a free 
streamline jet analysis, limited in its pertinence by the 
implicit geometric features of a jet flow originating in 
the horizontal (a = 180°) direction with no nozzle wall 
guidance on the ambient side. 

The most significant improvement has been afforded 
by the work of Strand, who generalizes the analysis to a 
free streamline model, restricted in generality and 
physical realism only by the introduction of an arti- 
ficial boundary condition requiring parallel flow normal 
to a line at the nozzle exit plane. This latter analysis 
is the first that obtains the realism, and offers quali- 
tative estimates of the effect of ground proximity in 
regulating the quantity of nozzle flow. Relevance is 
thereby maintained, without restriction, to zero alti- 
tude. 

A more complete analysis of the jet, including the 
effects of geometric detail of the nozzle, can be made 
to optimize the geometry of the nozzle and provide 
design data for estimating the actual operation of the 
system. In addition, it is of interest to study two 
classes of the problem. 

(1) Representative of local conditions in the curtain 
jet during climb or nonlevel flight, where some of the 
curtain jet air may flow into the support pressure re- 
gion—the bifurcated jet. 

(2) Representative of the steady-state level flight 
condition where, in general, there is no flow into the 
support pressure region—the deflected jet. 
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Fic. 1. A ground effects machine using a curtain jet. 


Case I—Bifurcated Curtain Jet 


We focus our attention on the detailed geometry and 
flow conditions at the two-dimensional nozzle as in 
Fig. 2a. The jet from the nozzle impinges on the 
ground surface and splits into two streams, one into 
the higher pressure support region and the other into 
the lower pressure ambient region. The boundaries 
of the efflux streams are typified by constant static 
pressure streamlines separating them from the still, 
constant pressure support and ambient regions. The 
configuration is defined by the inclination of the nozzle 
sidewalls (a), the altitude of the nozzle edge closest to 
the ground in comparison to the nozzle width (y,/f), 
and the relative alignment of the sidewall edges in ratio 
to the nozzle width 


5 = (v2 — 9;)/t (1) 


The analysis proceeds to evaluate the effectiveness 
of the jet in supporting a pressure differential and an 
estimate of the flow involved, as well as a detailed 
description of flow geometry, velocities, etc. 


Analysis of the Bifurcated Jet 


The situation has been reduced to a configuration 
amenable to the Helmholtz-Kirchhoff analysis sum- 
marized in standard hydrodynamic tests. The funda- 
mental entry for construction of the analysis lies in 
recognition of the fact that all flow boundaries are 
either (a) streamlines on straight solid walls, each 
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Fic. 2a. The bifurcated curtain jet configuration considered for 
analysis. 
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with uniform inclination, or (b) streamlines with con- 
stant velocity. The hodograph of the flow is there- 
fore bounded by (a) radial lines, and (b) circular arcs, 
Since the conjugate of the hodograph plane is a har- 
monic function, a set of conformal transformations 
(Fig. 2b) is possible to define the problem solution 
completely. 
The conjugate hodograph 
y=u-—iWw=qe”™ (2) 
is transformed into a polygonal area by the logarithm 
of its reciprocal : 


Q = In (U/v) = In (U/q) + 00 (3) 


This plane may, in turn, be transformed into the upper 
half-plane (¢) by the Schwarz-Christoffel transformation 
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Fic. 2b. The analytic mapping planes, bifurcated jet. 


d¢ 
aw=cf ~ : ———_ (4) 
(¢ — s)[(? — DE — ny(e — 7) ]}!”? 


In addition, the complex potential function (w) may 
be transformed into the same upper half-plane (¢) by 
Schwarz-Christoffel transformation. This gives, for 
the particular problem, 


(dw/d¢) = (Ut/r)(¢ — s)/(€2 — 1) (5) 


The location of the stagnation streamline in the half- 
plane is related to the relative strength of the two 
sinks by the ratios 


bo 


(U,4)/(Ut) = (1 — s) (6a) 


(Uot2)/(Ut) = (1 + s)/2 (6b) 


The mapping parameters (C, 7, 72, s) can be related to 
the physical geometry by appropriate integrations be- 
tween key points on the various boundaries. The 
local slope of the free streamline is related to the trans- 
formations by the relationships derived by integration 
along the free streamlines 


A(é) = —1Q(5) ka = —CG,(E) (7a) 
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ai In (U/U,) = Q(¢)J2 = CH, (12a) 
itn con- 

s there- In (U/U2) = Q(f)|™. = CH (12b) 
lar arcs, eile 
3 a har- 
mations df : 
. - support H, = | = ar : é >, (ida) 
solution Jo (¢ — s)[(f? — 1)(€ — n)(E — 7%) }!/? 
region 
® si dé 
F ambient H, = im = : - (13b) 
-o® (¢ — s)[(f? — 1)(€ — n)(E — 7) }!/? 
(2) region , 
Final definition of the geometric detail in terms of the 
‘arithm mapping parameters requires locating the nozzle edges 
in the physical plane (y; and y.). This can be done by 
C 
(3) c 
E 
- upper , i os , 
; Fic. 3a. The deflected curtain jet configuration considered for 
mation analvsis. 
0(€) = —iQ(¢) ii — x —CG(~) — r (7b) riba 
where 2 
cxccrrees fl . (p te di 3 
G(s) = ees —— ~ (da) 
J1 (¢ — s)[(E? — 1)(n — OE — 12) ]”? ~ 
“rr lo dé s 
smo 4’ G2(g) = Soe eT - (8b) 
J-1 (s — O(? — ID(n — OC — ve) ]}” 
The complete integration along the streamline to the 
nozzle tips gives 
A,(71) =-—-e = — CG,(r1) (9a) 
Ws 0 0.2 0.4 0.6 0.8 1.0 
6.(r72) = —a = —CGi(r) — @ (9b) 
= ? altitude ¥ 
. ‘ t 
t, from which we derive 
. . . Fic. 4a. Effect of nozzle angle on support pressure, bifurcated 
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= \ a bd . . - - . . 
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where 
lh = f sin 0,(¢ — s)dt/(¢2 — 1) (19a) 
1 
i= f sin 0.(¢ — s)d¢/(¢? — 1) (19b) 
-1 
Symmetric Bifurcated Jet 
The particular set of conditions s = Oand7 = —72, = 


r establish a symmetric picture in the mapping half- 
plane. This symmetry is reflected, then, in the 
physical plane and defines a perpendicular jet (a = 
90°) with equal static pressures, efflux velocities, jet 
widths, etc., on either side of the jet. Eqs. (8a), (8b), 
(13a), (18b), (19a), and (19b) are reduced readily to 
integrable forms. The problem is then a special case 
of the class of solutions treated by von Mises® and is 


summarized by the results 
U/U, = U/U2 = [(r — 1)/(r + 1)}!”? (20) 
t/t = te/t = [(r — 1)/(r + 1)]!/2/2 (21) 


yi/t = ye/t = (h/t) + 
[r/a(r + 1)] In [r + (7? — 1)1?] (22) 


Case II—Deflected Curtain Jet 


The bifurcated curtain jet geometry implies a net 
amount of flow into the high-pressure zone. In the 
physical situation of a ground effect machine in static 
equilibrium and where the pumping machinery does 
not draw air from the support pressure region, the net 
flow into the support region must be zero. A more 
pertinent model to this situation is the deflected cur- 
tain jet, Fig. 3a. In this case, both boundaries of the 
jet leaving the nozzle deflect to the lower pressure 
ambient region. The streamline on the pressure side 
joins the ground surface at a tangent point with no 
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Fic. 5a. Effect of nozzle angle on flow rate, bifurcated jet. 


stagnation region. Once again, the free boundaries 
of the jet are typified by constant static pressure 
streamlines separating them from the still, constant 
pressure support and ambient regions. The same 
geometric parameters (a), (y/t), and (6) define the 
picture, with the analysis aimed at evaluating the 
effectiveness of the jet in supporting a pressure differen- 
tial and estimating the flow involved and details of flow 


geometries, velocities, etc. 


Analysis of the Deflected Jet 


The case of the deflected curtain jet, as constructed 
in Fig. 3a, is also amenable to the Helmholtz-Kirchhoff 
type analysis. Referring to the set of transformation 
planes of Fig. 3b, the following analytic relationships 
may be derived. 

The (Q) plane, in which the flow field is represented 
as a rectangle, may be transformed into the upper half- 
plane (¢) by the Schwarz-Christoffel transformation 


O85) = CS a&/[(E -— PYG -— YD)? (23) 


The complex potential plane (w) can also be transformed 
into the (¢) plane by the function 
(dw/dé) = —(Ut/m)(a — 1)/(€ — a)(E — 1) (24) 
Continuity requires that 
Ut = Ug (25) 
As in the previous case, the mapping parameters 
(C, r, a) can be related to the physical geometry by 
integrations between key points on the various flow 
field boundaries. The local slope of the free streamline 
is given by 


A(é) = &(—&) = —iQ(¢)F = 
? r(g2 — 1 )1/2 (r? — 1)1/2 ; 
(—C/nisn| *. -3 | (26) 
E(r? — 1)}/2 r 


in closed form.? With the end condition, § = 7, & = 
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—a, we may solve for the constant 
C = (ar)/sn—[1, (r? — 1)'/?/r] (27) 
In addition, integration along the ground line gives the 


relationship between the two boundary streamline 


velocities 
Oe oy yt oF 
In (U2/Ui) = Q(6)]_, = (—2C/r)sn—"[1, (1/r)] (28) 


in closed form.’ Integration along one of the nozzle 
sidewalls relates the velocity upstream in the nozzle to 


the free streamline velocities 
ae aie 
In (U/U,) = Of) ], = 


; we = 
(—C/r)sn | =. <3 | (29) 


(a? — 1)*/? 
in closed form.’ Finally, the integration relating dis- 
tances in the physical (z) plane to the mapping half- 
plane implied by Eq. (15) may be used to locate the 
nozzle edges 


(y,/t) = (U/U)) [1 + (a — 1)h/r (30a) 


(y2/t) = (U/U2)(a — 1)Io/x (30b) 
where 

I, = f sin 6,d¢/(¢ — a)(¢ — 1) (31a) 

1 

7-—y 
Lh = | sin @.d¢/(¢ — a)(¢ — 1) (31b) 

e —1 

Results 


The results of the foregoing were evaluated nu- 
merically as detailed in the Appendix. 

The effect of the angle of inclination of the nozzle is 
presented in Figs. 4 and 5. The special case of a 
square-ended nozzle is used as a basis of comparison— 
that is, the situation in which a line connecting the ends 
of the nozzle sidewalls is perpendicular to the center- 
line. Specifically, we deal with the case 
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Fic. 5b. Effect of nozzle angle on flow rate, deflected jet. 
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6 = —cos a (32) 


Results are presented for nozzle angles (a) of 90°, 
112.5°, and 135°. We define the dimensionless param- 
eter 

(po ani pi) /(po sa Pi) =l1- (Us U;)? (33) 


to indicate the ratio of (the pressure differential avail- 
able for support of the ground effect machine) to (the 
dynamic head imparted to the curtain jet by the pro- 
pulsive machinery pumping nozzle flow from ambient 
conditions). 

Figs. 4a and 4b indicate that maximum support 
pressure at the higher altitudes is achieved by inclining 
the nozzle away from the vertical, pointing toward the 
support region. There is a diminishing rate of im- 
provement as 135° is approached, and for the bifurcated 
jet the effect is actually reversed, with maximum sup- 
port pressure effectiveness achieved at a nozzle angle 
less than 135°. 

Figs. 5a and 5b present data on the flow rate of the 
nozzles in terms of maximum rate of flow at the same 
supply and ambient pressures when the nozzle is far 
removed from the ground restriction: 


W/Wmar = U/U, (34) 


The data indicate that there is a consistent trend 
toward use of higher flows for larger nozzle angles for 
the bifurcated jet and the opposite trend for the deflected 
jet. The difference in trend is explicable in view of the 
increasing amounts of flow which are lost into the sup- 
port pressure region as nozzle angle is increased for the 
bifurcated jet, whereas no equivalent loss is possible in 
the deflected jet situation. 

Figs. 6 and 7 indicate the effect of relative alignment 
of the nozzle edges (6), using nozzles inclined at 112.5° 
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Fic. 6b. Effect of nozzle edge alignment on support pressure, 
deflected jet. 

as a basis of comparison. It is noted that retraction 
of the sideplate on the support pressure side of the nozzle 
results in substantial increases in available support 
pressure for the bifurcated jet (Fig. 6a) and in the 
opposite trend for the deflected jet (Fig. 6b). The 
nozzle wall, A-F (Fig. 3a) in the deflected jet sustains 
a positive pressure gradient in the flow direction. 
This suggests that flow separation may be experienced 
on the wall in the physical situation of a real fluid. 
Such premature separation would in effect reduce the 
effective length of the wall, so that Fig. 6b is, therefore, 
indicative of the trend toward decreasing pressure 
support effectiveness associated with separation on the 
nozzle wall. In all cases, retraction of the nozzle wall 
on the support pressure side results in some reduction 
of flow (Fig. 7), although the effect is not significant 
enough to plot in the case of the deflected jet. 

The most significant conclusion derived from com- 
parison of the two types of jets is the profound effect 
which extraction of fluid from the pressure side of the 
jet has in deterioration of curtain jet performance. 
Typically, for a support pressure function of 0.70, the 
bifurcated jet rises to only 13 percent of the altitude 
of the deflected jet apparatus, while using 48 percent 
of the flow. 

It should be noted that, although the analysis is 
two-dimensional, it is readily referable to the annular 
jet for configurations in which the width of the jet is 
small compared to the diameter of the annulus. In 
addition, it also seems reasonable to use various values 
of the solution to represent /ocal jet conditions for non- 
axisymmetric cases of noncircular machines, and ma- 
chines in motion with flight stagnation pressure on 
the forward side, and machines in nonlevel flight. 

Various engineering parameters for optimization 
and analysis of particular configurations, involving lift 
associated with jet reaction as well as support pressure, 
are readily deducible from the presented data. 

Some experimental treatment of the curtain jet is 
available in the literature. Von Glahn* has investigated 
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an annular jet equivalent to the analytic case (a = 
90°), over a range of values (y,/f), for three values of 
6. Unfortunately, no data were taken in the range 
where the altitude of the nozzle from the ground line 
is of the same order as width of the jet (y, < ¢), which is 
probably the most pertinent regime of the analysis. 
Von Glahn’s results do point out the limitation of a 
two-dimensional analysis by establishing the regime 
for larger altitudes where the annular jet coalesces 
and causes significant entrainment and under-pressure 
in the base region of the nozzle. 


Appendix 
Numerical Evaluation 


Except for the special case of the symmetric jet 
cited in the text, the integrands of Eqs. (Sa), (Sb), 
(18a), and (13b) are not integrable in closed form in 
terms of simple functions. A numerical solution was 
carried out, employing the following devices. 

Referring to a typical integral, Eq. (Sa), the linear 
variable transformation 


m = (¢ — 1)/(n — 1) (Al) 


was made to convert integration range to 0 < m < 1. 
With this transformation, the integral has the form 


m 
Gi(m) = f F(m)dm/([(1 — m)m]}'/? (A2) 
0 
In order to cope with the improper behavior of the 
integrand at m = 0 and m = 1, actual numerical 
integration was performed on the equivalent form 
G,(m) {- F(m)dm F(O)dm F(1)dm 
ni(m ee oe 
o [(1 — m)m|}}/? m}!2 (1 — m)}!? 


+ 2F(0)m'/? — 2F(1)[(1 — m)”? — 1] = (A3) 


( Continued on page 871) 
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Supersonic Shear Flow Past an Airfoil 
Between Two Parallel Walls 


INGE L. RYHMING* 


Convair, A Division of General Dynamics Corporation 


Summary 


The supersonic flow with assigned Mach number gradient in 
the span direction past a straight wing between two parallel 
walls is studied using the small-disturbance theory. The govern- 
ing equation for the disturbance pressure on the airfoil, together 
with the boundary conditions on the airfoil and at the walls, is 
solved by the method of separation of variables. Upon sepa- 
ration the problem is reduced to a Sturm-Liouville eigenvalue 
problem and to the solution of the telegraph equation. 

As an application, a certain Mach number profile is selected 
and the resulting pressure distribution on a parabolic arc airfoil 


is computed. 


Symbols 
Sy: 9: Cartesian coordinates 
l = wing span 
h(x) = wing thickness 
p(x, y, z) = disturbance pressure 
U(y) = velocity of oncoming flow 
a,(y) = speed of sound of oncoming flow 
pi = pressure of oncoming flow 
oily) = density of oncoming flow 
Mi(y) = Mach number of oncoming flow 
w = z-component of disturbance velocity vector 
Y aa Cp/Cy 
z = 2(M,? — 1)1/2 
F(x, Z) = function satisfying Eq. (9) 
(yy) = function satisfying Eq. (10) 
d = separation constant 
S(A) = spectrum function 
T(x) = Bessel function of order n 


Laplace transform variable 


Gn, bn, Cm = constants 


5(x) = Dirac delta function 

I'(n) = gamma function 

g(é) = kernel function defined by Eq. (21) 
t, n(t) = variables defined by Eq. (26) 

bu = constant defined by Eq. (28) 

T = constant 

Mo = Mi(y = 0) 

A, BC constants 


(1) Introduction 


| gpiuamencnd Honda! presented the solution to the 
problem of incompressible shear flow past an air- 
foil between two parallel walls. The supersonic counter- 
part to this problem is of interest in compressor and 
turbine design as well as in wind tunnels. The solu- 
tion to this latter problem can, as will be shown, be 


easily obtained. To solve the problem, the general 
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method of attack will be essentially the same as that 


used by Honda in his paper. 


(2) Governing Equations and Boundary 
Conditions 


The small-disturbance equation for the flow with 
assigned velocity and density gradients past a thin 
airfoil is given in reference 2. In the present case, the 
coordinate system is chosen as shown in Fig. 1. The 
x-axis is in the streamwise direction and the y-axis is 
along the leading edge of the airfoil. The two parallel 
walls are given by y = 0 and y = /. The airfoil is 
situated between x = 0 and x = 1. Gradients in ve- 
locity and density are assumed only in the y-direction. 
The small-disturbance equation for the flow past the 
airfoil can then be written as 


(M,2 — 1)prz — M,?(M~*d,)y — Pez = 0 (1) 


where M, = M,(y) is the Mach number of the oncoming 
flow and p is the disturbance pressure. The Mach 
number distribution along the y-axis is determined by 
M, = U(y)/a\(y), where U(y) is the velocity and a;(y) 
the speed of sound of the oncoming flow a,*(y) = 
¥ [P:/ p(y) | 
Introducing the transformation 

2(M,? — 1)” = 2 (2) 

in Eq. (1), we have: 


(M,? _ 1) (Pez =— i = M{?(M, "Dos = V0 (3) 


Let z = h(x) describe the wing surface; then the 
boundary condition on the wing surface can be written 
as 

(w,)e-90 = (d*h/dx*)U(y) for0O< x < 1O<y</l 


(4) 


where w is the z-component of the disturbance velocity 
vector. The z-component of the linearized momentum 
equation is 


U(y)w, = —[1/pily) |p. (5) 
Eq. (5), together with Eq. (4), gives for the boundary 
condition : 


(pz)zm0 = —U?(d*h/dx*?) = —ypiM,*(d*h/dx*) 
With the transformation given by Eq. (2), this is 


(Pz)se0 = — VP: Mi?(M,? — 1)~-'? d*h/dx? (6) 
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M,(y)>1 








Fic. 1. Coordinate system for an inviscid supersonic shear flow 
past an airfoil between two parallel walls. 


The boundary condition which is to be fulfilled at the 
walls is 


by = 90 for y=0,/ (7) 


Also Pb = p,2 = pb: = 0 for 


(3) Solution 


To solve Eq. (3) with the boundary conditions Eqs. 
(6), (7), and (8), we apply the method of separation of 


variables. Assume that an elementary solution of 


Eq. (3) can be written as 

p(x, y, 2) = F(x, 2) Wy) 
Then we get [Eqs. (9) and (10)]: 
F,, = Fy + X*F = 0 (9) 
d/dy((1/M,")(dy/dy)| + (Mi? — 1)/MYy = 0 (10) 
where \ is the separation constant. The general solu- 


tion is formed as 


p = > S(A) F(x, 2; ) Wy; A) (11) 
r 


where S(A) is the spectrum of eigenvalues. Eq. (11), 
together with Eq. (6), now gives 
y~'p,-1M,—? (M,? — 1)? > S(A) [Fs(x, 2; X) ]s-0 X 
¥(y; A) = — dh/dx? (12) 
In view of the boundary conditions we have 
F3(x, 2; \)z-0 = —(d*h/dx*?) on 0 Sx < 1 (18) 
F= F,=F,=0 for x — jz)}<0O (14) 
y—1p1 1M, -? (M,? — 1)*/? x SA)v(y; A) = 1 for 
O<y<l (15) 
d/dy[¥(y;\)] = 0 ony =0,/ (16) 
With these boundary conditions, the solution of Eqs. 


(9) and (10) is uniquely determined. 


(a) The Function F(x, z; \) 
The problem of integrating Eq. (9) together with 
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the boundary conditions Eqs. (13) and (14), has been 
given by Chang and Werner.* Their solution of the 
general form of the telegraph equation is obtained by 
the Riemann integration method. With the actual 
boundary conditions, the Riemann method does not 
give an explicit solution formula for the function F, 
Instead it gives an integral relationship between 
F(x, 2) at any field point and F;;(x, 0) on the boundary. 
The problem is thus reduced to the solution of an 
integral equation, which in our notation reads 


—h"(x) = f g(é) Jo[A(x — &) |dé (17) 
0 


In reference 3, g(&) is obtained by Laplace transforma- 
tions, and the result is shown to be valid for \ < 1. 
In our case, however, a solution for \ > 1 is needed. 
For \ > 1 it is found convenient to modify the result of 
reference 3 slightly; thus, taking the Laplace transform 
of Eq. (17), we have 


Lig) = L(—h")(s? + d2)"? (18) 


To cancel the branch points introduced by the factor 
(s? + \?*)!/? in the inverse transform of g, it is assumed 
that the function h” may be represented by the following 
series expansion: 


ro) 


—h"(x) = > [anx"J, (Ax) + b,x" +1J,,(Ax)] + 
n=0 
ul 


> Cmi(x — tm) (19) 


m=0 


where 6(x — ¢,,) represents the Dirac delta function. 
The last term is included to provide for a finite number 
of ordinary discontinuities of # at points x = t,. The 
coefficients c,, are determined by the magnitude of these 
discontinuities. To determine the coefficients a, and 
b, we first use the multiplication theorem‘ of the Bessel 
functions 

1 


v! 


Tni(X x) = DS [1/2 x (1 — A)” Jug (x) (20) 
) 


y=( 
Then expanding the J,;,(x) in power series of x and 
comparing term by term with the Maclaurin expansion 
of h”(x), the series Eq. (19) can be found. Introducing 
Eq. (19) in Eq. (18) and carrying out the inverse trans- 
formation, the following expression results for g(é): 
wv 
' Dn 
g(t) = a 6(é) + Ls —— X 
n=1 I'(n) 
fanV(n + 1/2)E"-¥? Tne (A E) + 


V2 Jn-12 (A &)] + 


i) | (21) 


The solution of the boundary-value problem can now 
be obtained, and with the analysis in reference 3 the 


V2V | 
[1/n b, T(n + 3/2)& 


n=0 I'(m) 
M 
Ji{vX(E ad tm) | 
po Cm r 7 7 + 6 (€ — 
m=0 é —_ Ton 
—E>tm 


solution is 


—F(x,0;’) = 1/A { g(&) sin [A(x — &)]dé (22) 
0 
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(b) The Function y(y; » ) 

Eq. (10) for the function ¥(y; A), with the boundary 
conditions Eq. (16), represents a Sturm-Liouville 
eigenvalue problem. There are infinitely many eigen- 
values and the functions ¥(y; A) form an orthogonal set 
in (0, 7). With the notation ¥(y; A,) = y, the condi- 
tion of orthogonality is 


Ser (23) 


1 
f (M2 — 1)/M\*v,~,dy = 0, 
0 
To normalize we set 
l 
Vt [ (M,? — 1)/M,? yp,” dy l (24) 
0 


The spectrum function, S(A,), is determined from 
Eq. (15). Multiplying Eq. (15) by ¥(y; An) vpi X 
(M,? — 1)'/*/~', integrating both sides over y, and 
applying Eqs. (23) and (24), we get 


l 
Sn) = vp if W(V; An) (Mi? — 1)/*2 dy (25) 
0 


(4) Example 


In order to elucidate the theory we shall now in- 
vestigate a simple flow model. A special Mach number 
distribution J/,(y) is selected by introducing the trans- 


formation 
t= [ [M,?(v) — 1]? dy, Wly; A) > n(t; A) (26) 
0 
in Eq. (10). We get 


dM, ; M,? — 92 
dy M,(M,2 — 1)? 


" 


, 9 = 
n +Ay=0 = (27) 
We now choose the particular case of Mach number 
profiles characterized by 


dM, a 2 _ 
dy MM, — 1)?” 


= 2u = constant (28) 


Eq. (27) then takes the simple form 


n” — 2Qun’ + r\*» = 0 (29) 


A 
M(y) 


Fic. 2. The function 2 tan! (M,? — 1)!/2 + (M,? — 1)73/2 = 
0.0726 y + 2.5 


SHEAR FLOW 863 





Cp./M6-! 
~ 8h 
. t-y=3 
lk 
y=0 
TWO-DIMENSIONAL y =O 
0} >, 
-| 
| 
Fic. 3. Local pressure coefficients cp = 2p/ypiMi%(y) as a 
function of x at the two walls y = 0, 3. 


The Mach number distribution is obtained by integra- 
tion of Eq. (28) and is 
2 tan—!(M,? — 1)!/2 + (M,? — 1)7? 

2uy + const. (30) 
This Mach number distribution is also used in references 
2and 3. A graph of this Mach number profile is shown 
in Fig. 2. It is seen that M, = +~/2 divides the profile 
in two different intervals of practical interest. The 
eigenvalue problem can now be solved. The boundary 
conditions to Eq. (29) are 


n’(t;4) = 0 ont=0,7 (31) 


where 7 is the value of ¢ for y l. Ford < pitis seen 
from Eqs. (29) and (31) that the only solutions are 


const. ) 


n=C, X9=0 (C5 


(32 
n=0, \=u ” 


To determine the constant C, and S(A = 0), we must 
transform Eqs. (24) and (25) to the ¢ representation. 
This can be done if we know the variation of Mach 
number with?. Indeed we can obtain this variation by 
combining Eqs. (26) and (28), with the result 

Qut = In[M,2(My? — 1)'/2/My?(M,? — 1)1/?] (33) 
where My is Mi(y = 0). With Eqs. (26) and (33), 
Eqs. (24) and (25) are 


I-1My-2(Mg? — 1)2! i ne" dt = 1 (34) 
/70 


S(A) = yp,/l | n dt (35) 
JV 
The constant C and .S(0) are now determined as 
C? = QulM,?(M,? — 1)-/2 (1 — s"™")-1 (36) 


(Continued on page 884) 








Plastic Buckling of Longitudinally Stiffened 
Plates' 


H. L. SUJATA* 


Atomics International** 


Summary 


Formulas are developed which define the plastic buckling 
coefficients for a simply supported plate stiffened by longitudinal 
stiffeners and subject to uniform compression in the direction 
of the stiffeners. The equations are developed by consideration 
of the strain energy of plastic buckling. 
in detail are the plate stiffened by one, two, and an infinite 
It is shown that design charts of the 


The cases considered 


number of stiffeners. 
buckling coefficients can be readily developed for very long 
stiffened plates. An approximate method is presented which 
permits evaluation of the plastic buckling coefficient for stiffened 


plates with any aspect ratio. 


Symbols 
a = length of the plate in the X direction 
b = distance between stiffeners 
¢ = aye + 7 
d = 32 + Y 
e = B/E, — 1 
f = ase + 7 
g = ans t+ 7 
h = aus + Y 
l = integer indicating the number of one-half sine waves 


in the X direction 
j - integer indicating the number of one-half sine waves 
in the Y direction 


k = buckling coefficient = o,,(b*t)/(m?Dx ) 
a minimal value of the buckling coefficient 
ki = kB?/(1 — v?) 

m = integer 

n = integer 

p = a3 + 

q = 143 + Y 

r = integer defining the stiffener location 

t = thickness of the plate 

u = total number of bays 


w = deflection in the Z direction 


A = plasticity coefficient 

Ai; = a particular Fourier coefficient 
Amn = Fourier coefficient 

As area of a stiffener 

B = plasticity coefficient 

B, = rigidity of the stiffener = E,/, 
D = plasticity coefficient 

Dx = Ft?/(12(1 — v?)] 

E = ela8tic modulus of elasticity 
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oF = tangent modulus of elasticity 

E, = secant modulus of elasticity 

F = plasticity coefficient 

Ip = moment of inertia of a plate unit width = (°/12 
K = ratio of the plastic k»;, to the elastic Rmin 
N = stiffness of the plate = E/, 

P = load acting on a stiffener Aor 

R = total number of stiffeners 

Vi = internal energy 

Ve = external energy 

Vis = internal energy of a stiffener 

Ves = external energy of the stiffener 

Vip = internal energy of the plate 

Vep = external energy of the plate 

Qiju = Ait+ 2(B + 2F)(1j6/u)? + D(78/u)4 
Bx = ratio of applied stress = o,/o, 

B = aspect ratio = a/b 

Bmin = value of 8 corresponding to k»,i, 

6 = ratio of areas A,/bt 

7 = ratio of stiffness = B,/bN 

Yeq = equivalent elastic stiffness ratio 

v = Poisson's ratio 

E A 4 
Or, Oy = unit stresses 

Cor = critical unit stress 

x = A+ 8(B + 2F)6 2+ 16Dp6' + 7 
y, = terms in the equations for the plasticity coefficients 


Introduction 


7 PROBLEM Of elastic buckling of stiffened plates 
has been studied by several authors.'~> Timo- 
shenko appears to have been the first to investigate 
this problem when he used the strain energy approach. 
Barbré® investigated this problem by solving the 
differential equation of plate buckling with the ap- 
propriate boundary conditions incorporating the effect 
of the stiffeners. Bijlaard? recently applied his 
method of split rigidities to the problem. The most 
comprehensive study, particularly from the design 
point of view, is the work of Seide and Stein’ which 
develops charts for plate-stiffener design. Cox and 
Riddell® have made an excellent analytical presentation 
of this problem. 

The behavior of stiffened plates in the inelastic 
region was analyzed by Schleicher.’ Recently Koll- 
brunner and Meister’ have presented Schleicher’s 
method as an approximate solution to this problem. 
One of Schleicher’s important assumptions was that 
the plate remained isotropic in the inelastic range and 
that only a modified plate stiffness value was needed 
to include inelastic behavior of the plate. Bijlaard’ 
has shown that the behavior of a plate during inelastic 
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Fic. 1. Idealized stiffened plate. 


buckling is in better agreement with experimental 
results if the plate is considered anisotropic. 

The following discussion uses the inelastic stress- 
strain relationships developed by Bijlaard’” to de- 
termine the buckling coefficients of plates with longi- 
tudinal stiffeners. The results may be incorporated 
into charts which can be used in design procedures. 


General Remarks 


The structural component which is analyzed is 
shown idealized in Fig. 1. The following assumptions 
are incorporated into the subsequent work. 

(1) The area and moment of inertia of the stiffener 
are concentrated at a point—that is, at the intersection 
of the center line of the stiffener with the middle plane 
of the plate. While this assumption leads to a more 
simplified analysis, this structural element is generally 
constructed with the stiffener on one side. In this 
case, it is necessary to determine an effective width of 
the plate which will add to the moment of inertia of the 
stiffener. !!_ An analysis of the plate effective width 
in the inelastic range has been made and is presented 
in reference 12. 

(2) The stiffener has no torsional rigidity. 

(3) The stiffener does not buckle locally before 
buckling of the plate-stiffener combination can occur. 

(4) The shear deformation and flexibility of the 
stiffener are neglected. 

(5) The stiffener and plate are loaded with a uniform 
compressive stress. 

The buckling load of the stiffened plate shown in 
Fig. | is determined by making use of the strain energy 
concept. This concept states that buckling will occur 
when the internal energy of the system equals the work 
done by the exterior load when the system is given a 
small lateral displacement. Stated mathematically, 


Vi = Vz 
or more specifically, 
Vie + Vis = Ver + Ves (1) 


The expressions for the external and internal energies 
of the plate have been given by Bijlaard.* These 
can be written as 
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Vep 


to a bu ow 9 
= f [ —) dxdy (2) 
m4 0 J0 Ox , 


_ { dw \?}* 
1F ( y] dxdy (3) 
Oxdy , j 


The expressions for the internal and external energies 


of the stiffener are: 
Ed, (° [o*w\? 
(RY we 
2 0 Ox? y=rb 


oA * [ow\? 

r ae 8 ~ 

Ves =— f ( dx (5) 
2 0 Ox y=rb 


It will be noted that the expression for V;s uses the 
tangent modulus of elasticity instead of the elastic 
modulus. This is consistent with the results of 
Shanley,’* who showed that column buckling occurs 
at the Engesser load. 

The deflection surface of the plate is given by the 
double sine series in the X and Y directions: 


Vis 


and 


w= >), Dd Ann sin (max/a) sin (nry/ub) (6) 
1 


m=l1n= 


This expression satisfies the boundary conditions 


that the bending moment and the deflection on all 
edges be zero. 
The deflection of the stiffener is given by Eq. (6) 


but with the stiffener location substituted for y. This 
gives: 
w= > DY Ann sin (max/a) sin (nar/u) (7) 


When Egs. (6) and (7) are substituted into Eqs. (2), 


(3), (4), and (5) and the resulting expressions substituted 
into Eq. (1), the following equation is obtained: 


* The coefficients A, B, D, and F are functions of the applied 
stresses and the ratio E,/E and E,/E and are defined by: 






A= W/o B= be/tn D = vs/ty F = 1/(2 + 2h + 3e) 
where 
vi = (1 — 264)?E + 3(1 + e)E, 
yo = (2 — Be)(1 — 26e)E — [(2 — Ba)(1 — 28x) — 
4un? — 3eBx| Ey 
v3 = (2 — Be)?E + 3(1 + e)Be2E, 
vs = [(5 — 4v)(1 + Be?) + 2(5v — 4)Be + 3en2]E — 


(1 — 2v)[(1 — 2v)n? — 3(1 + e)Bs JE; 
in which B+ = o,/o, = ratio of applied stress 


9 = Be* — Be + 1 
e = E/E, - 1 


The above relationships are derived in reference 10. 
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Nrj oes m*n?B? Dn 4B x . . nar\? 
er J Uu > pe Ag,* | Am‘ + 2(B + 2F) —— + 4- 2>> Yr ,? m4 p>» Amn Sin l 
b°B? i m=ln=1 Uu- u4 y=) m= n=] Uu f 
tor = ; oe Se a —— (8 
ast s _ nrr\? 
Uu L > Amn?m? + 2 } » 6, » m?* (x Ay, sa ) 
m=l1n=1 r=1 m=1 n=1 u 
The summation on r = 1 to r = R sums into the : ; : ; 
; : Pe ; symmetrically about the stiffener. Symmetrical buck- 
equation the number of stiffeners. The resulting : ; ae ; 4 os 
i ti ers is bie ling will occur if the stiffener does not possess sufficient 
expression is therefore in its most general form. ; : ‘ 
‘ aan stiffness (E,J;) to keep the plate from deflecting F 
For any given case, all of the terms of Eq. (8) are ; : : 
iy ee ie outward at the point at which the two are connected. é 
known with the exception of A,,,, the Fourier coefficient. : : ; ; c 
ae : Antisymmetrical or nodal buckling will occur if the i 
In order that these coefficients be such that they yield ht age i ‘ae F 
-—F ; “ty stiffener is sufficiently rigid to prevent the movement. c 
the minimum buckling load, the 0¢,/0A,; is taken Sage : ; : ‘ 
; aa This implies that the stiffener acts as a knife-edge 2 
and the resulting expression set equal to zero. ' 5 
, "Mage : - support and that it allows rotation to take place : 
When this is done, the following expression results: ‘ i : : a 
freely when buckling occurs. If the plate is quite long : 
R -) e : . o 7, 
on ee and nodal buckling takes place, there is a fourfol ¢ 
UA 430434 + 2 : y7rt* sin (jar/u) z. x : ; g : I . : - 3 
; r= ea increase in the critical stress when this structural 3 
kp? element is compared with an unstiffened plate with 
A in Sin (nar/u) = ~~ ee tes I yt* + the same overall dimensions. 
R ce For the plate which is reinforced but where E,J/, is 
2>> 6,1? sin (jrr/u) +> A in Sin (nar | (9) not sufficiently large to cause nodal buckling, it will 
r=1 =] . ° ° ° 
* be shown that there is still an increased capacity over 
where k = (o¢,b*t)/(m?Dx) the unstiffened plate. 
and aij, = Ait + 2(B + 2F)(ij8/u)? + D(jB/u)* ; ; 
Symmetrical Buckling 
Eq. (9) is a linear equation in general terms which ; ; 
be (9) é a : “8 By using Eq. (9) with R = 1 and u = 2 and con- F 
will be used in the evaluation of the buckling coefficients weed! d ‘ ; ; , 
; . ’ 2 sidering only the first term of the series defining the 
for a plate with one, two, and an infinite number of . ; : anes ect 
5 deflected shape in both the X and Y directions, the 
stiffeners. ; ; : oe ‘ 
following expression for the critical stress is obtained 
when 7 = j = 1 are substituted: 
Plate With One Stiffener ‘ , 
Ocr = ka*Dx/b*t (10) 
In this case two buckling modes are possible. The : 5 
it? gigs where Ds = EI,/(1 — ») 
plate may buckle symmetrically, resulting in the 
displacement of the stiffener, or it may buckle anti- and of 
5.0 ——— r I “ 
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Fic. 2. Buckling coefficient vs. aspect ratio. Plate reinforced by one stiffener. 
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Fic. 3. Compressive buckling coefficients for a simply sup- 
ported long plate with one longitudinal stiffener. 


| l 
- (B+ 2F)p? + 6 DB + ¥ 
) 


A+ 


(1 + 6)B" 
(11) 


The minimal value of k and the corresponding value 


of B»in can be found as: 


j l(c + d)(1 + 6) — 2y6| — V 


Nodal Buckling 


When the stiffener is sufficiently rigid to prevent 
outward movement of the plate at the center, the plate 
will buckle about the stiffener. In this case, the plate 
buckles into two one-half sine waves in the Y direction; 


therefore, Eq. (9) is used with ¢ = 1, 7 = 2. This 
yields 
Ocr = ka*Dx/b*t 
where 
k = (1—»){[A + 2(B + 2F)6? + Dp‘)/e2} (15) 


The buckling coefficient for nodal buckling is inde- 
pendent of the number of panels in the element when 
using the notation of this paper. Consequently, the 
nodal buckling coefficient is the same for all of the 


subsequent cases. 


Symmetrical and Nodal Buckling Coefficients 


The curves for the buckling coefficient k vs. the 


‘ 
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l 
2(A + ¥) +5 (B+ 2F)Bmin? 


Renta _ (1 aid py?) 4 = = 
(1 + 0)Bmin* 


where 


Buin” = 16(A +- Y) D 


defines the value of 8,;, corresponding to the first 
minimum point of the buckling curves. 

The expression for k given as Eq. (11) yields results 
which are generally sufficiently close for engineering 
purposes when 8 = Bmin. However, when 6B < Byin, 
more accurate results are obtained when the first two 
odd terms of the sine series in the X direction are con- 
sidered. When this is done, the following two simul- 
taneous linear equations result: 


kB? 
c— (1+ 4) . Ayn - 
1 — p* 


(1 tore 
Ay = 0 
(1 — p*) 


where 


l 
A + 5(B + 2F)B? + (Dp*/16) +-4 


9 
A+-= 


S] 
d= (B + 2F)p? + 16 Dg* + ¥ 
) 
Setting the determinant of the coefficients A,, and 
Aj3 to zero and solving for k, the following expression 
is obtained: 
(c + d)(1 + 6) — 2y6]? — 4(1 + 26)(cd — y?)) 


(14) 
2(1 + 26)B? f 


aspect ratio for one particular set of plate-stiffener 
combinations and for the various states of plasticity 
are shown in Fig. 2. It will be noted that the curves 
representing the plastic state are of the same form as 
the curve for the elastic case. 

It should be pointed out that y is defined as E,/,/ bEI, 
and that this stiffness ratio is constant for the entire 


range of values of E,/E. This means that as E,/E 
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Fic. 5. Compressive buckling coefficients for a simply supported plate with two longitudinal stiffeners. 
decreases, J; increases proportionally, giving the the solid curves are for symmetrical buckling. Addi- 


constant value for y. It may also be thought of by 
regarding the stiffener as elastic for all values of F,. 
The plate becomes weaker as -;// and E,/E decrease, 
but the rigidity of the stiffener remains constant. 


Buckling Coefficients for Stiffened Plates of Infinite Length 
With One Stiffener 


The presentation of curves for k vs. 8 which could 
be used for all possible combinations of E/E, E,/E, 
y, and 6 as well as over the commonly used range of 
values of 8 would be extremely voluminous. Instead 
of this, the minimal values of k for symmetrical buck- 
ling for all possible values of -,/E and E,/E were 
computed. These values of k could be conservatively 
used for all values of 8 but they would be very con- 
servative when 8 < £,,;, on thez = 1 curve. A method 
which gives reasonably good results in this range is 
presented in the Appendix. 

It was found that the values of k for symmetrical 
buckling varied only slightly with E,/E.'*  Conse- 
quently, if one were to use the minimum value of k 
for a given set of values of E,/E, E,/E, y, and 6, the 
resulting buckling coefficient would be conservative 
and little error would result therefrom. 

The values of & for nodal buckling, however, are 
direct functions of both £,/E and E,/E; a simplifica- 
tion similar to that for symmetrical buckling is not 
possible without extremely conservative results. 

The minimal values of k for nodal and symmetrical 
buckling are given in Fig. 3 for the case where 6 = 0.2. 
The dashed curve corresponds to nodal buckling and 


tional cases are considered in reference 12. 

With computed values of y and /,//, one can employ 
Fig. 3 to determine the buckling coefficient. For 
example, if the computed value of E,/E = 0.7 and y = 
40, it is seen that the plate will buckle nodally if 
E,/E = O or 0.1, and symmetrically for all other 


values. 
Plate With Two Stiffeners 


The plate which is reinforced by two stiffeners can 
be analyzed in a manner similar to the foregoing 
analysis of a plate with one stiffener. 

In this case, the stiffeners are connected to the 
plate at the one-third points and Eq. (9) is modified 
for this case by substituting RK = 2 and u = 3. 

Three buckling modes must be investigated in order 
to establish the coefficient k. These three modes 
symmetrical, antisymmetrical, and nodal—are shown 
in Fig. 4. 


Symmetrical Buckling 

As an initial approximation, the plate is assumed to 
buckle into one one-half sine wave in the Y direction. 
By using only the first term in the series defining the 
deflection, there is obtained: 


2 D 
A+, (B+ 2s + — B+ 4 
= — yp? 6) 
k (1 v*) (1 5)B2 (16 


When the first two terms in the series are considered, 
the following expression for & results: 
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fl(g + A)(1 + 6) — 276] — V I(g + h)(1 + 6) — 2y6]? — 4(1 + 26)(gh — y?)) 


k = (1 — »?): - 
\ b 2(1 + 26)B? ( (17) 


[he minimal value for k can be readily found from 

































































































Eq. (16) as , . : a . i 

. buckling mode is possible. The buckling coefficient 

Rmin = [2(4 + vy) + 9 (B + 2F)Bmin?|/[(1 + 8)Bmin?] corresponding to symmetrical buckling generally will 

: be considerably lower than that corresponding to 

(18) antisymmetrical buckling. However, to establish that 

where Bmin = BV(A + y)/D (19) the plate will not buckle antisymmetrically for all 
values of 8 required further investigation. 

As in the case of the plate with one stiffener, the Barbré* found that. antisymmetrical buckling is 
value of k given by Eq. (16) yields satisfactory results possible for elastic plates with two longitudinal stiffen- 
when 8 2 8B,,;, of the curve 7 = 1. However, when ers when 6 is small. However, Seide and Stein® 
8 < Bmin, the results are unconservative, giving a reported an error in Barbré’s results and concluded 
coefficient higher than the actual value. In this range, that antisymmetrical buckling was not possible in 
Eq. (17) can be used to give good results. the elastic case. In the plastic case, however, where 

: ‘ ' the plate behavior is governed by different constants, 
Antisymmetrical Buckling at ; : : . cal ‘ 

it is desirable to investigate the possibility further. 

The question arises as to whether an antisymmetrical The equation for the antisymmetrical buckling 
coefficient can be found as: 
r e e / ©19 ' 9 

b lb + gC + 6) — 276] — V [(6 + 9) + 46) — 2y6)? — 4(1 + 286)(pg — vy’) 7" 

= (1 — p?): (20) 
2(1 + 26)8? 

Eq. (20) is plotted in Fig. 5 for the case where y = 2, 
$§= 04, E/E = 09, E,/E = 0.8 and £,/E = 0.3, is unlikely to occur with the small ratios of stiffness. 
E/E = 0.1. Two other cases were investigated!*—that is, where 

It will be noted that for the computed cases, buckling y = 100 and with the same ratio of moduli. In these 
is never antisymmetrical. Inasmuch as the cases cases, it was found that the buckling coefficient corre- 
investigated are near extreme values for the ratio of sponding to antisymmetrical buckling was greater 
the moduli of elasticities E,/E and E,/E as well as than that corresponding to nodal buckling. 
near the lower extreme for the ratio of stiffness y, it is The explanation of this phenomenon is that in these 
reasonable to conclude that antisymmetrical buckling cases the plate has little rigidity relative to the stiffener 

5.0 ] 
’ | 
ELASTIC NODAL | 
3 BUCKLING | 
4.0 E | 
™ — PLASTIC NODAL BUCKLING ELASTIC SYMMETRICAL Ey; I< | 
= Eox= Ey = BUCKLING awe 
s Ff %209 Ey " Y= pep 720 | 
= “ p 
2 3.0 | 
ud “ As 
oO . PLASTIC S= br 700 
ts. = SYMMETRICAL BUCKLING 
tn ™ Es/=0.9 —t,-=0.8 | 
o = 
© 20 
2 . 
ain + 
x - PLASTIC NODAL BUCKLING 
Se > | Es/e =0.3 ' Ee =0.! 
a) C eS PN 
LO + 
. PLASTIC i=l 
sl SYMMETRICAL BUCKLING 
s Ee =0.3 —Y¥e=0.! 
oo Cttbit tt non ee Pe poi tii sicaisiealiadiaanal 
0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 


ASPECT RATIO (2 = B) 


Fic. 6. Compressive buckling coefficients for a simply supported plate with an infinite number of longitudinal stiffeners. 
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and the plate is readily forced into nodal buckling. 
However, to force the plate into antisymmetrical 
buckling would mean that the stiffeners would have to 
be displaced and this would require a higher stress 
and increased energy. It may be concluded that as 
the stiffener becomes more rigid relative to the plate, 
the less likely antisymmetrical buckling becomes. 

Still another way of reaching the same conclusion 
is to base the explanation on the conclusion of Seide 
and Stein® for the elastic case. It was pointed out 
previously that the behavior is as if the stiffener re- 
mains elastic for all values of -,/E and E,/E. 

However, the resistance of the plate diminishes with 
decreasing moduli ratios and the ratio of rigidities is 
the smallest when the plate and stiffener are elastic. 
Consequently, if the plate does not buckle antisym- 
metrically in the elastic case, it will not in the plastic 
case—as defined herein. 


Design Charts 


The minimal values of k for nodal and symmetrical 


k <= (1 sis vy 


Explanation of Buckling Curves 

Fig. 6 shows a typical set of curves for a plate 
reinforced by an infinite number of stiffeners. These 
curves were computed by Eq. (22). 

It will be noted that there is little difference to the 
buckling coefficients even when the plate is well into 
the plastic range. The reason for this is that the 
effect of the supports parallel to the applied stress is 
very small. Consequently, the buckling stress of a 
panel is directly proportional to the composite radius 
of gyration of the stiffener and the plate about the 
center of the plate. 

It will be noted from the definition of the ratio of 
stiffness E,/,/bET, that the curves of Fig. 6 are plotted 
This states that the rigidity of 
However, as the plate 
This is 
becomes 


for y being constant. 
the stiffener /,/, is constant. 
becomes more plastic it becomes weaker. 
analogous to postulating that the plate 
thinner while remaining elastic. 

Since the buckling stress of the plate-stiffener com- 
posite element is directly proportional to its radius 
of gyration, it is reasonable to expect the buckling 
stress to change very little as the plate becomes thinner 
because it is at the center of the composite section. 
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buckling can be incorporated into a series of curves, 
as was done for the case of one stiffener. 


Plate With an Infinite Number of Stiffeners 


Two buckling modes are investigated for this plate 
stiffener combination. The first case is that of sym- 
metrical buckling with the entire plate buckling in 
one one-half sine wave in the Y direction. The other 
case is that of nodal buckling. 


Symmetrical Buckling 

As an initial approximation to the buckling surface, 
it will be assumed that the deflected shape is defined 
by the first term of the infinite series. 

It can be shown that the buckling coefficient for 
this case becomes: 

k= (1—»){(A + y)/[ + 587]! (21) 

A better approximation to the buckled surface is 
made by using the first two applicable terms in the 
series defining the plate deflection. 

The following equation for the symmetrical buckling 
coefficient is obtained after some manipulations: 


flé + x). + 8) — 2y8] — VUE + XC + 8) — 2yd)? — 4(1 + 2)(EX — YY (22) 
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Appendix 


A Procedure for Evaluating the Buckling 
Coefficient for all Values of the Aspect Ratio 


A family of curves similar to Fig. 3 can be obtained 
for any case desired by using Eqs. (12) a dn(15) for 
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PLASTIC 


one stiffener and Eqs. (18) and (15) for two stiffeners. 
These curves give the correct value of the buckling 
coefficient only if the aspect ratio 8 is large. This 
buckling coefficient can be used, nevertheless, for all 
values of 6 and the results will always be conservative. 

However, it will be noted from Figs. 2 and 5 that for 
the smaller values of 8, a great deal of difference exists 
between the coefficients for nodal buckling and sym- 
If 8 were such that the point fell 
and 


metrical buckling. 
on the buckling curves (Rk vs. 8) where 7 = 1 
8 < Bmim, the design buckling coefficient might be 
quite conservative. 

In order to extend the usefulness of the method 
contained herein, and to eliminate this ultraconserva- 
tism, a method of evaluating k in this range is given. 
The procedure which follows is applicable to the cases 
where there are one or two stiffeners. 

This particular method notes and makes use of the 
similarity in the curves of Seide and Stein for the 
elastic range and those given in Figs. 2 and 5. It 
consists of the following procedures: 

(1) Compute an equivalent elastic stiffness ratio 
Yeq by 


ve = [y(1 — v)]/(E./E) = (0.91y)/(E./E) (23) 


2) Use y¥-q¢ and 6 to locate the proper elastic curve 


in reference 5. The value of 6 (ratio of areas) remains 
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from the elastic curve corresponding to y,, and 6 the 
minimum value of the buckling coefficient. If the 
value cannot be determined directly from the curves, 
it can be readily computed by: 


“/ 
Bmin uv | + Yeq (24) 


and 
Rmin - 


(3) Select the plastic value from the pertinent chart 
(Fig. 3) and divide it by the elastic ki, found in 
paragraph 2 above to determine the multiplicative 


[2(1 + Yee) + (Bmin?/2)]/[(1 + 8)Bmin?] (25) 


constant 


K = Plastic ky» j,/Elastic Ry» i (26) 


(4) The nature of the design problem will prescribe 
8B. Using this value of 8, read off the buckling coeffi- 
cient k,, from the curves of reference 5 which correspond 
Multiply this value by the constant K 
The resulting value 


tO Yeq and 6. 
computed in paragraph 3 above. 
is the approximate buckling coefficient which can be 
used in the subsequent design calculations. There 
will be some instances where y,, will be larger than the 
largest elastic y plotted in reference 5. It 
though, that the values of y,, falling into this group 
are not very large and the above method covers most 
of the cases where there is a large difference in the 
buckling coefficients of nodal and symmetrical buckling 


appears, 
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(Continued from page 860) 
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Che Curtain Jet 

which converts the integrand to finite form at m = 0 


and m = 1, permitting numerical treatment. 

Numerical integration was carried out using the 
trapezoidal rule with Am = 0.02. 

Integration of Eqs. (19a), (19b), (81a), and (31b) was 
also carried out in numerical form, using the results of 
the integration of Eqs. (7a), (7b), and (26) to evaluate 
the inclination (@, and #2) of the free streamlines. 

The integrand of Eq. (19a) is properly behaved at its 
upper limits, but is indeterminate at the lower limit, 
in the form (0/0). Using the linear transformation 
of the independent variable as in Eq. (Al) above, a 
first-order form of the integral Eq. (19a) for small 
values of m is represented by 


Lhlo ~ —C[2m/(1 — r)] + (A4) 


This solution was used in the first half of the first incre- 
mental integration interval, that is, 0 < m < Am/2, 
instead of the trapezoidal rule. An analogous expression 
was developed for the integration of Eqs. (19b) and 
(3la) in the region of the lower integration limit. No 
special problem of indeterminateness was encountered 
in numerical integration of Eq. (31b). 

The numerical results using the above procedures 
were accurate approximations to the closed-form 
solution for the symmetric bifurcated jet to within 
+(0.001 over the entire range of variables studied. 


Solution was carried out for a systematic variation 
of the variables 7, 72, and s (bifurcated jet) and the 
variables 7, a, and a (deflected jet) over a wide range 
of their possible values. Solutions for dependent 
variables were interpolated from the results at values 
of independent variables of engineering interest. 
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Direct Calculation of Pressure 
Distribution on Blunt Hypersonic 


Nose Shapes With Sharp Corners' 


MAURICE HOLT* 
University of California, Berkeley 


Summary 


The method of Belotserkovskii!~* for calculating hypersonic 
flow fields past a circular cylinder is extended to deal with axially- 
symmetric flow past sharp-cornered nose shapes, in particular, 
flat-headed cylinders. Results on 
In the present paper Belotserkovskii’s 


spherical segments and 
spheres are also included. 
first approximation is considered, and comparison of calculated 


pressure distribution and shock shape with experimental results 


shows very good agreement. 


Symbols 

r,0 = polar coordinates 

M = Mach number 

q, P, p = dimensionless total speed, pressure, and density 
(Velocities are made dimensionless in terms of the 
maximum velocity, pressure and density in terms 
of stagnation pressure and density respectively. 
Other symbols are defined in the text.) 

u,v = velocity components in polar coordinates 

7 = adiabatic index 

k = (7 — 1)/2) 

¢ = p/p? 

€ = distance from body to shock measured along radius 
vector 

y = stream function 

x = angle between normal to shock and axis of symmetry 

T =(1 — g?) V4 = 

t = Tv 

h = Tu 

3 = puv 

g = 2kp + pv? 

H = kp + pu? 

Indices 


0—body surface 
1—behind shock 
*—sonic point 
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(1) Introduction 


7” DIRECT METHOD for calculating the flow field 
around a blunt-nosed body in a hypersonic stream, 
due to Belotserkovskii, has been used and discussed 
extensively.1~> In the existing formulation of the 
method it is assumed that the body contour is smooth, 
and it is difficult to treat bodies with discontinuities in 
slope. 

In the present paper the method is extended to apply, 
in the first approximation, to hypersonic flow past 
nose shapes terminating in sharp corners—in particu- 
lar, the forward face of a flat-headed cylinder and 
spherical caps. 

The spherical caps are to be’attached to some after- 
body and their geometry is restricted so that sonic 
conditions will be attained at the sharp shoulder. This 
will be the case provided that the angular extent of the 
cap is less than that of subsonic flow on a smooth sphere 
(about 40°) and also provided that the discontinuity 
of slope at the joint of the cap to the afterbody is not 
small. This discontinuity must further be large enough 
to ensure that the flow field past the nose is independent 
of the flow field near the afterbody. The precise condi- 
tion for this is given in reference 6. These restrictions 
are always fulfilled on flat-headed cylinders. 

As in the case of spheres or other axially symmetric 
bodies, the motion is determined by three first-order 
ordinary differential equations for the shock detach- 
ment distance, the shock angle, and the surface ve- 
locity. The new application, however, requires a 
change in the downstream boundary conditions. On 
the body, in place of the condition that derivatives of 
physical variables remain bounded through the sonic 
point, using this property to locate its position, we now 
require that sonic conditions are attained at a fixed 
point on the body and allow the derivatives to become 
infinite. Off the body, derivatives across the sonic 
line must be bounded, as in the case of smooth con- 
tours. 

Calculations are made on a spherical cap of semi- 
angle 26.21° at incident Mach numbers of 2.01 and 
4.76. The results are compared with corresponding 
experimental data’ and excellent agreement is found. 
For spherical caps the equations of motion are most 
conveniently expressed in terms of polar coordinates, 
based on the center of the sphere. The same equations, 
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combined with the condition of smooth transition 
through the sonic point on the body, determine the flow 
field round a sphere. For completeness, therefore, re- 
sults are also presented for the sphere, at Mach num- 
bers of 5.8 and 18. At the lower value, these compare 
well with Oliver’s experimental data.® 

To treat the flat-headed cylinder Belotserkovskii’s 
first approximation is expressed in terms of Cartesian 
coordinates (with the x axis along the front face). 
Apart from this change the calculation is similar to that 
for spherical caps. Complete details are given by Gold 
and Holt. The results, obtained at a Mach number 
of 5.8, agree well with experimental data obtained by 
Boison and Curtiss.’ 

In contrast to results for smooth bodies, the applica- 
tion of Belotserkovskii’s method to nose shapes with 
sharp corners determines the flow field only up to the 
radius or coordinate line through the sonic corner. 
On certain smooth bodies, in particular those for which 
the body curvature is nonincreasing in the downstream 
direction—e.g., spheres, smooth sphere-cone and sphere- 
cylinder combinations, and prolate  ellipsoids—the 
method can be applied successfully a considerable dis- 
tance downstream of the sonic point. On bodies with 
increasing curvature, such as oblate ellipsoids, the 
method can only be carried a short distance beyond the 
sonic point. But since the flow field is supersonic here 
its computation can be continued by the method of 
characteristics. 

To continue the flow field calculation beyond the 
sharp sonic corner on the present shapes it is necessary 
to join the Belotserkovskii solution to a solution in 
series in the centered expansion round the corner. In 
the case of Belotserkovskii’s first approximation the 
series solution given in reference 11 is sufficiently 
accurate. If a higher order approximation is applied 
to the present shapes (a difficult task still remaining 
to be attacked), the solutions should be matched to the 
detailed and exact series solution derived by Vaglio- 
Laurin.'* Such series solutions carry the flow field 
around the corner and provide initial data on a char- 
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Fic. 1. Polar coordinate system. 





Shock shapes for a sphere. 


Fic. 2. 


acteristic line terminating the expansion adjacent to 
the afterbody. The flow field can then be continued 
along the afterbody by the method of characteristics. 
This aspect of the calculation depends on the shape of 
the afterbody and is not considered in the present 


paper. 


(2) Equations of Motion in Polar Coordinates 


We consider the motion of a perfect gas with con- 
stant specific heats past a sphere or spherical segment. 
Far upstream the motion is uniform and supersonic, 
and there is a detached shock wave ahead of the body. 
Following Belotserkovskii we introduce dimensionless 
variables, referring the velocity to the maximum veloc- 
ity, the pressure and density to the stagnation pressure 
and density ahead of the shock, and lengths to the radius 
We use the system of 


of sphere or spherical nose. 
in which 7 is the 


polar coordinates shown in Fig. 1, 
distance from the center of the sphere and @ is the angle 
between the radius and the axis of symmetry. 

The motion can be determined from the momentum 
equation in the 7 direction, the continuity equation, 
the equation of conservation of entropy, and Bernoulli’s 
equation. The first three of these may be written 
as follows. 


Momentum: 


(0/dr)(r?H7 sin 6) + 
(0/00)(rs sin 8) — grsin@=0 (1) 


Continuity: 
(0/dr)(r?h sin @) + (0/06)(rt sin 0) = 0 (2) 


Entropy: 
et 2% = p/p” = $V) (3) 
where y is the stream function and 
H = kp + pu? Ss = puv | 
g=2kp+ pw? h= ru t= 7 > (4) 


r=(1—g)YO-> k= (y— 1)/29) 
Bernoulli’s equation gives 


p=(1-—@)p (5) 
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Fic. 3. Pressure distribution on a sphere. 


and hence, from the entropy equation, 
= ia = 2) ; 
p = To (6) 


The boundary conditions are that the body surface is 
one of constant entropy, and hence of constant ¢, and 
that the shock equations are satisfied across the de- 
tached shock wave. 

We introduce as dependent variables the distance « 
between the body and the shock, measured along the 
radius, and the angle x between the tangent to the shock 
and the vertical. These are shown in Fig. 1. 

We now integrate Eqs. (1) and (2) by Belotserkov- 
skii’s method. Using the first approximation we as- 
sume that the dependent variables s, g and ¢ vary 
linearly with 

£ = (r — 1)/e 
between the body and the shock. The problem is 


then reduced to the integration of the following three 
ordinary differential equations for ¢, x, and vp: 


(de/d6) — (1 + e) tan (0 — x) = 0 (7) 
a(dx/d0) + b(de/d@) +c = 0 (8) 
l — |] 
sin 6 ce + 3e) (7 = (2 - a) x 
y—W\y +1 
‘ / Lv dx de 
i--_° rr 7 ~+d "= 0 9 
; ) a’ *n**ar? (9) 


where a, b, c, d, e, and f are the following known func- 
tions of 6, €, x, and vw: 


a = sin 0 (2e? + 3¢e) O(pityv)/Ox 


b = —sin 0 (Ze + 3) pir 
c = 6(1 + ©)? sin A(kp; + pitts?) —6Rpo sin 06 + 
ra) v 
cos 6 (2e? + 3e) pyyvi sin 6 (2e2 + 36) - aoe _ 


sin 0 (e? + 3e)(2kpo + poo”) + (Qe? + 3e) X 
(2kpi + piv”) 

d = sin @ (2e? + 3¢e) O(717;)/O0x 

e = sin @ [(2e + 3)vo7. — (2e + 3)r1%4] 

f = 6(1 + ©)*7,.m sin 0 + (2e? + 3e)712 cos 6 + 


: ae 4 O(7121) - 
sin 0 (2e? + 36e) oe + cos 0 (e? + 3¢) Tov 
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The index 0 denotes values on the body. The values 
of quantities behind the shock, denoted by index 1, 
are given in terms of x by the shock equations. De 
tails of these are given by Holt‘ and also by Traugott.5 

At this point an important difference between the 
present and Traugott’s® paper should be noted. Here 
the dependent variables s and ¢ are themselves assumed 
to vary linearly with 7, and this appears to be con- 
sistent with Belotserkovskii’s original approach. In 
Traugott’s treatment this assumption is made about 
the weighted variables vs and rt. This difference 
has a noticeable but not large effect on results for 
spheres. The effect could be expected to be more 
pronounced in the second approximation, especially in 
the case of nonspherical shapes. 


(3) Boundary Conditions and Integration of the 
Basic Equations 


Eqs. (7), (8), and (9) must be integrated along the 
shock layer between the axis of symmetry and the 
neighborhood of the sonic line with boundary condi- 
tions applied at both ends. The boundary conditions 
at the axis of symmetry are the same in all problems, 
but those near the sonic line depend on whether the 
body is a complete sphere (Case 1) or a spherical seg- 
ment (Case 2). 

Case 1 

In this case the flow along the body is free from 
singularities, and we require that all flow variables shall 
be continuous in @ with finite derivatives. In Eq. (9) 
the coefficient of dvj/d@ vanishes when vp attains the 
sonic speed v* = (y — 1)/(y + ) ae Hence, dv,/dé 
can remain finite at this point, @ = @*, only if the re- 
maining terms in Eq. (9) vanish there. Accordingly, 
the boundary conditions in Case 1 are as follows. 


Ate’ = 0: 


) 
eEe= @& f (10) 

At @ = 6*: 
K = d(dx/d@) + e(de/dé) + f = 0 (11) 


The shock detachment distance ¢ is, of course, not 
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Fic. 5. Pressure distributions for a spherical cap. 


known in advance and must be so chosen that condition 
(11) is satisfied. 


Case 2 

In this case the body has a sharp corner at the junc- 
tion of the nose and the afterbody. At this corner 
sonic speed will be attained, and the flow will turn the 
corner through a centered expansion, so that dvo/dé will 
be infinite at the sonic point. The nature of the singu- 
larity is determined from Eq. (9), which states that 


(v* —T Vo) d(v* — V%), d@ 


remains finite and nonzero as @ approaches 6* (the 
asterisk refers to sonic conditions, as before). It fol- 
lows that, near the sonic point, 


(v* =— Vo) ~~ (@* , 9)? 


Similar arguments apply in the case of the flat-headed 
cylinder. This approach to sonic velocity in proportion 
to the square root of distance from the sonic point is 
dictated by the assumption that ¢ and s vary linearly 
across the shock layer. It is a close approximation 
to the exact law of approach established by Guderley,!® 
in which the square-root law is replaced by a two- 
fifths power law. In a higher order Belotserkovskii 
approximation it may be possible to introduce Guder- 
ley’s condition by treating a small neighborhood of the 
sonic corner separately from the rest of the field, to 
which Belotserkovskii’s method is applied. This re- 
finement is not expected, in itself, to make a great deal 
of difference in the numerical results. 
The boundary conditions are now the following: 


At 0 = 0: 
sie a 
At 0 = 6: 


where @, is the known value of @ at the shoulder of the 
segment. 
In this case, less subtlety is required in finding the 


PRESSURE DISTRIBUTION 875 


correct value of ¢. Sonic conditions will be attained 
upstream of 6 = 0, if the estimate of «is toolow. The 
only difficulty that arises is in using a difference process 
near the singularity in du/d@. To overcome this we 
stop the normal difference process for integrating Eq. 
(9) as soon as dv/d@ becomes too large, and extrap- 
olate up to the value v = v* by means of the formula 

v* — uw = A(o* — 0)” (14) 
determining A and @* from the earlier part of the 
integration. The correct value of « corresponds to 
6* = 6, in Eq. (14). 

The numerical work was carried out, in all cases, 
on an IBM 650 computer, which is slow by present 
standards. The running time for sharp-cornered 
shapes was shorter than the corresponding time for a 
smooth contour. In the former case a converged 
solution, giving all unknowns to within five-figure 
accuracy, could be obtained in ten iterations. On the 
IBM 709 computer the calculation of one case, for a 
sphere or sharp-cornered nose shape, can be completed 
in less than five minutes. 


(4) Results and Comparison with Experiment 


The method is applied in Case 1 to calculate the 
flow past a sphere in a uniform stream of Mach num- 
bers 5.8 and 18. At the lower Mach number results 
are compared with experimental results derived by 
Oliver. At Mach number 5.8 the specific heat ratio 
is taken to be 1.4. At Mach number 18 the value 1.2 
is used. 

The results are shown in Figs. 2 and 38. Fig. 2 shows 
the shape of detached shock wave, and Fig. 3 the sur- 
face pressure distribution. It is seen from Fig. 2 that 
the calculated shock shape at \J = 5.8 agrees very 
closely with the shape observed by Oliver.2 The calcu- 
lated pressure distribution is below the Newtonian 
curve (which is to be expected) but still reasonably 
close to Oliver’s experimental points. More calcula- 
tions on smooth spheres or sphere-cone combinations 
are given by Traugott.® 

In Case 2 calculations were made on a spherical cap 
of semi-angle 26.21° at Mach numbers of 2.01 and 4.76. 
Figs. 4 and 5 show shock shape and pressure distribu- 
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Fic. 6. Shock shape and pressure distribution for a flat-headed 
cylinder. 
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tion at these two Mach numbers. Results are com- 
pared with corresponding experimental results derived 
at the Jet Propulsion Laboratory.’ Agreement is 
excellent and shows that Belotserkovskii’s method is 
reliable, even in the first approximation, on sharp- 
cornered bodies. The experimental results confirm 
the existence of a singularity in the velocity gradient 
at the sonic point. 

Fig. 6 shows the shock shape and pressure distribu- 
tion for a flat-headed cylinder, calculated at a Mach 
number of 5.8. These curves are compared with experi- 
mental data obtained by Boison and Curtiss,!° extrap- 
olated from different Mach numbers to the theoretical 
value. The only measured shock detachment distance 
is on the axis. 
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A Theory of the Two-Dimensional Laminar 
Boundary Layer Over a Curved Surface’ 


K. T. YEN* ano K. TOBA** 
Rensselaer Polytechnic Institute 


Summary 


The purpose of this paper is to present a theory to account for 
surface curvature effects on the two-dimensional boundary-layer 
flow which approaches a potential flow at free stream. 

The problem of two-dimensional viscous flow is first formulated 
by using the streamlines and their orthogonal trajectories as the 
generalized coordinates. A boundary-layer approximation is 
applied to the Navier-Stokes equations and the Gauss equation 
in the generalized coordinates to yield the boundary-layer equa- 
similar solutions of the 
By a simple 


tions. The conditions under which 
boundary-layer equations exist are determined. 
transformation, the governing differential equation can be ex- 
pressed in a form which reduces to the Falkner-Skan equation 
for zero surface curvature. 

Numerical results for a similar solution which corresponds to a 
flow over a curved surface with zero surface pressure gradient 
have been obtained. The velocity profiles in the boundary layer 
and the wall skin-friction distribution for concave and convex 
surfaces are presented. The wall skin friction for a convex wall 
is found to be higher than the Blasius value for a flat plate. On 
the other hand, for a concave wall, the skin friction will drop be- 
low the Blasius value as the curvature increases, but it appears to 
reach a minimum, and beyond this minimum point it will increase 
The same flow problem was treated by Murphy? by a 
Comparison of Murphy’s results 


again. 
different method of analysis. 
with those obtained by the present method reveals some basic 
differences in the boundary-layer characteristics. In particular, 
Murphy’s results indicate that the wall skin friction for a convex 
surface is smaller than the Blasius value, while for a concave wall 


it is higher. 


Symbols 

C; = skin friction coefficient 

c = a measure of surface curvature 

he, h,~! = metric measures along é- and 7-directions, respec- 
tively 

heo = value of hg at » = 0 

H = defined by Eq. (31) 

R = a typical Reynolds number 

R:, Ry radii of curvature of é- and 7-directions, respec- 
tively 

u,v = velocity components in x- and y-directions, 
respectively 

UE, Uy = velocity components in §- and 7-directions, re- 


spectively 
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a = Cartesian.coordinates 
ay) = 28-1 
7B = wedge angle 
f wre” 
»* = R'/% 
6 velocity direction 
A = hy/hz 
& 7 = generalized coordinates 
y¥ = stream function 
w = vorticity 
ig 
x am rv ldg 


Introduction 


ers BOUNDARY-LAYER THEORY is generally 
assumed to be applicable to flows over curved 
surfaces if the surface curvature is small,' but no 
surface-curvature effect is actually considered. In 
reference 2, Murphy obtained the differential equations 
and boundary conditions for a boundary-layer flow over 
a curved surface with either ‘“‘moderate’”’ or “‘large’’ 
curvature. He applied his theory to the flow over a 
particular surface with zero pressure gradient to deter- 
mine the velocity profiles and the wall skin-friction 
distribution for various surface curvatures. 

Although Murphy’s analysis appears to be generally 
reasonable, the accuracy of his numerical results may 
be questioned. The main difficulty may be due to the 
method of analysis adopted by Murphy, which led to 
the problem of determining a function satisfying a 
third-order, nonlinear, ordinary differential equation 
and the boundary conditions f = f’ = 0, » = 0, and 


f’ > 2/(1 + 2An), 9 > @ (A is a surface-curvature 


parameter). The third boundary condition is evi- 
dently a difficult one to handle. The differential equa- 
tions, on the other hand, have the artificial singularity 
at Ky = —1 introduced by the method of analysis. 

In reference 3, a new boundary-layer theory which 
takes into account the surface-curvature effects was 
proposed by the first author of this paper. First, the 
Navier-Stokes equations are transformed into general- 
ized coordinates, which are chosen to be the streamlines 
and their orthogonal trajectories. The boundary-layer 
approximation applied to the Navier-Stokes equations 
and the Gauss equation in the generalized coordinates 
follows Prandtl’s argument concerning the orders of 
magnitude of the terms in these equations (see also 
reference 4). Two boundary-layer theories were pre- 
sented. The first one (P-approximation) is essentially 
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Fic. 1. Transformation of coordinates. 


equivalent to Prandtl’s boundary-layer theory. The 
second one (C-approximation) takes surface curvature 
into account. 

The conditions under which similar solutions of the 
boundary-layer equations in the generalized coordinates 
exist have been obtained in reference 3. It is found 
that for the C-approximation, the surface curvature 
and the surface velocity distribution of the potential 
flow are restricted in a manner analogous to that of the 
wedge flows in the conventional boundary-layer theory. 
Two ordinary differential equations and appropriate 
boundary conditions have been obtained. 

Numerical results for the similar solution, which 
corresponds to a flow over a curved surface with zero 
surface-pressure gradient (the same as Murphy’s), have 
been obtained by using IBM 650 computers. Details 
of the numerical analysis can be found from reference 5. 
The purpose of this paper is to present the theory de- 
veloped in reference 3 and the numerical results ob- 
tained in reference 5. 


Basic Relations in Generalized Coordinates 
Let &(x, y) and n(x, y) be two functions with 
(O&/Ox)(On/Ox) + (OF/Oy)(On/Oy) = 0 (1) 


so that the two families of curves § = constant and yn = 
constant are orthogonal to each other. Let h;~'dé and 
h,—'dyn be the elements of length along the &-lines and 
the n-lines, respectively. One can show that® 


ox\?2 dy\2]-1/2 

7 ts F Nas 
ai Ox oy); [4 ” 
is (3s) . (55) | . 


oh; a Oh,, ec , 

h, on = —h,R; ’ h; dt = —~AR.- (3) 

Ss 

where &; and R, are the radii of curvature of the é-lines 

and the 7-lines, respectively. The two functions /; and 
h, satisfy the Gauss equation 
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0 (h, Oh; 0 /h; Oh, 
at, ag hie ? =—1 = 0 (4) 
On Az” On O& a" Of 7 


The Navier-Stokes equations for two-dimensional 


incompressible flows in the generalized coordinates ¢ 


and 7 can be written in the following form :’ 


Oo7l . Ow 
h; or 9 v~ no p —- v.08 = (h, IX) On (5) 


j O i - (h./R Ow (G) 
1 Oy -* +p} — yw = (A;/K) >: ») 


where v? = »,? + v,”, and 
fo. a) ) 
» = hh,4\— (v,/h,) — v;¢/h) 7 
W . 8 lde n ” on (  . { 

The continuity equation is 
(0/0E) (v;/h,) + (0/0n)(v,/h;) = 0 (7) 


In the above equations, v; and v, are the velocity com- 
ponents along the &-lines and the 7-lines, respectively, 
p is the static pressure, and w is the vorticity. All these 
quantities are dimensionless. is a Reynolds number 
based on a typical velocity, a typical length, and the 
kinematic viscosity. All the above geometrical and 
flow quantities are taken as continuous and differenti- 
able whenever necessary. 

If w and v are the velocity components in the Car- 
tesian coordinates x and y, 


u = v,h,(Ox/OE) + v,h,(Ox/On) (8) 
v = v,h(Oy/0E) + v,h,(Oy/On) (9) 
When potential flows are considered (w = QO), or- 


thogonal curvilinear coordinates can be chosen such that 
E(x, y) and n(x, y) are the real and imaginary parts of an 
analytical function. Thus, 


h,=h,=h (10) 


Eq. (4) shows that In / is a harmonic function. 
Let (x, y) be the stream function such that 


dy = h,(O¥/On), %y= —he(Oy/dE) (11) 


and the continuity equation is automatically satisfied. 
Let the streamlines be used as the n-lines, and 7 = y. 
It follows that z; is the velocity magnitude, v, = 0, and, 
from Eq. (11),% = #,. Introducing the function 


A = h,/h; = v;/hy (12) 
and eliminating p from Eqs. (5) and (6) yields 
36 (ae) ty (ay) Bae 
where 
w = —h,*d(Or/On) (14) 


The Gauss equation, Eq. (4), can be written as 


re) (; 0 In *) ‘ fe) (x° In **) 0? (5) (15) 
ae fs — {\— = —|- 5 
Of \A OE On On Of? \A 


For a potential flow, \ = 1, Eqs. (13) and (15) are 
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satisfied automatically. Thus, if a viscous flow over a 
rigid surface, on which y is chosen to be zero, approaches 
a potential flow at free stream, the proper boundary 
conditions to be used for A are 


\=Oaty=0 (16) 
A>lasy~>o (17) 
and, in view of Eq. (14), 
OA\/On > Oas Po = (18) 
The boundary conditions for h; may be chosen to be 
h; = known value at y = 0 (19) 
and, from Eq. (3) 
\(Oh,/0n) = —R;-' at y = 0 (20) 


Eqs. (12) through (20) (with, possibly, some appro- 
priate initial conditions) can be considered as a mathe- 
matical formulation of the two-dimensional viscous- 
flow problem. This formulation evidently constitutes 
an inverse approach, but is found to be convenient 
for analyzing the surface-curvature effects of boundary- 
layer flows. A direct approach would involve a given 
body, and would require correction of the outer, 
potential flow for displacement-thickness effects and 
other comparable phenomena, as discussed briefly in 
the Concluding Remarks of this paper. 


Boundary-Layer Approximations 


The boundary-layer approximation is carried out in 
the generalized coordinates by introducing the trans- 
formation 


q° = RN% (21) 


From Eqs. (13) and (15), the first-order boundary-layer 
equations are obtained 


do /. , on d a2 /, ,ar\ o 
h;? = r h,? (22) 
O& \ * On* On* L On* On* 


(0/On*) [A(O In h,/On*)] = 0 (23) 


The next boundary-layer equations are of the order 
R-. The boundary condition (20) becomes 


NOh,/dn*) = —(R¥?R)-, n* = 0 (24) 


When R; at n* = 0 is of order 1, the boundary condition 
(24) should be applied to the first-order boundary-layer 
equations. Thus, the surface curvature effects would 
appear in the first-order boundary-layer approximation 
(the C-approximation). 

When the radius of surface curvature, R; at n* = 0, 
is of the order RK’? or larger, the surface-curvature 
effects will be, at most, of the order R~'. Thus, within 
this approximation, it is permissible to take h; as a 
function of £ only, and Eq. (23) is satisfied automati- 
cally. Eg. (22), the only equation needed to be con- 
sidered, becomes 





| . d f. da? 4 
h,*(&) = h,*(&) r (25) 
dE On* eM On* \" On*? 


and the boundary conditions are 


r 0, »* 0) (26) 
A> 1, n*—> (27) 
Od/On* > 0, *—> @ (28) 


The above approximation can be called the P-approxi- 
mation, because it is essentially the same as Prandtl’s 
boundary-layer theory as shown in reference 3. 

Now use of similar solutions for the C-approximation 
is considered. Let the similar parameter be 


¢ = n*/g(é) (29) 


where g(£) is to be determined. The solutions of Eqs. 
(22) and (23) are assumed to have the form 

ACE, n) = ACS) (30) 

h,(é, n) h(—)H(¢) (31) 


Substituting Eqs. (30) and (31) into Eq. (22) yields 


AE d»? _@ (ie) d E d (we) | 
= de 7 dt dt dt dé dé 


(32) 
where 
ay = (g*/hyo”) (d/dE) (hyo*/g) (33) 
and 
a, = g(dg/dé) (34) 


A necessary condition for the existence of the similar 
solutions is that a, and a». should be constants. Thus, 
it is found, with az taken to be 1, 

g = (2&)"/? (35) 


and 
hg = gure” (36) 


Substituting Eqs. (30) and (31) into Eq. (23) yields, 
after an integration, 


a | Surface - YC 











ee “= “ 0 2] 2 

Y| = - 

5 | OG 

+) Cory ol | 
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f 4 } + 2ST ———— a J 
| == +} 02 | 
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Fic. 2. Surface shapes, a; = —1. 
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\(dH/dt) = CH (37) 


where C is a constant of integration. From the bound- 
ary condition (24) and Eq. (37), the surface curvature 
is found to be 
] 1/2 1/4 
’ (ai — 1) /4 2Q 
(R;-!),*-0 = —C ( R) & (38) 


9 


Hence C is a measure of the surface curvature. 
From the above analysis, the differential equations 
governing the similar solutions for the C-approximation 


are 
dy? d , ar? d d , dr? 

aT” — ¢ de (we ~) = de as. (we 7) | (39) 
\(dH/dt) = CH (40) 

The boundary conditions are 
1=0, ¢=0 (41) 
H=1, ¢=0 (42) 
A> 1, [> o@ (43) 
drx/di >0, (> o@ (44) 


The similarity variable is 
f = 9*/(2g)¥? (45) 
When the similar solutions exist, it is seen that 


v% = hy, = hp(E)A(O)AS) (46) 


, ; 
where 7 = exp {cf (dg¢/X) t If the potential flow 
0 y 


prevails everywhere, then A = 1, and /j,o(£) is seen to be 
the surface velocity. Therefore, similar solutions for 
the C-approximation exist when the potential fiow 
yields a surface-velocity distribution of the form ¢'*"/4 
and the surface-curvature distribution is given by 
Eq. (38). Existence of potential flows satisfying these 
conditions is assumed here. 
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Fic. 3. Velocity profiles. 
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It is of interest to note that the differential Eq. 
(39) can be transformed to a convenient form which 
reduces to the Falkner-Skan equation’ for zero surface 
curvature. This derivation is carried out in the 
Appendix. 


Surface Shapes Yielding Similar Solutions 


The shape of the surfaces which yield similar solutions 
can be determined as shown in reference 5. From 
Eq. (38) and the following geometrical relation® 


h,(00/d&) = —1/R; (47) 


where @ is the angle between a é-line and the x-axis 
(Fig. 1), it is found that along the surface (ny = 0) 


6 = % + C(2Ré)'/? (48) 
The surface coordinates can be obtained from 
Ox/O0E = (1/h;) cos 0, Ov/OE = (1/h;) sin @ (49) 
and 
we wt 
— -(a+)/<, E _ -—(ai+1)/4 -: - 
x = | é cos édé, y= | E sin 6dé 
0 0 
(50) 


If the distance along the body surface is denoted by s, 


aed - 
dé 4 ee . 
s -| - 2 aie its 
0 heo 3 — ay 


Along the 7-lines, the distance 7 is given by 


"dy 2 aan [ 7H dy | 
_ = eo A exp 4 —-¢ -~> dt 
’ J h,r R 0 bite o XS ” 
9 1/2 od 9 ike 
= +(1—ai) 4) ’ as | a 
= E t J l — ex — ( 4 f ye 
(aa) aes | J, r |} = 


Thus, the similarity variable may be taken as (2/2)'/” 
x ns (lan /(3 an) | 

The constant % in Eq. (48) can be determined as 
follows. When C is zero, the function HT becomes 
The above solutions should reduce to 


then 


identically 1. 
those of the wedge flows in the Prandtl theory.’ If the 
wedge angle is 78, the potential-flow velocity will be 
of the form s”, m = 8/(2 — 8). By using Eqs. (36) 
and (51), it is found that 


a +1 = 28 (53) 

and it is reasonable to choose 
6 = (1/2)B8 (54) 
The surface shape for the case a; = —1 is shown in 
Fig. 2. Since the surface velocity 4, becomes constant 
for a, = —1, this case corresponds to a flow without 


surface-pressure gradient. The surface curvature is 
given by 

1/R; = —C(R/2é)/? (55) 
Since s = from Eq. (51), this flow is exactly the same 
one studied by Murphy.? The surface-curvature pa- 
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Fic. 4. Boundary layer velocity profiles near the wall. 


rameter C is equal to —2'/?A, where A is the surface- 
curvature parameter used by Murphy. When C = 0, 
the problem reduces to the boundary layer flow on a 
flat plate. 


Numerical Example 


The problem is to determine the function \, for which 
a single differential equation can be obtained from 
Eqs. (39) and (40): 


d?2(d*n2/de*) + [ACA + AE + (1/2) (dd*/dE) | d?A?2/dg? + 
(4C? + 2C¢ — aigf)dd?/dé = 0 (56) 


The boundary conditions are Eqs. (41), (43), and (44), 
which suggest a trial-and-error method of solution by 
starting from the wall, ¢ = 0. 

The velocity profile v; vs. n(R/s)'/? can then be 
found by using Eqs. (46) and (52). The wall skin 


friction is given by 


9 1/2/dy2 8 
C,(Rs)'/? = ( — ) ( ) (57) 
3 Qa) d¢ r=0 


where C, = 2 u(Ou/Oy)y0/pU* and (U = hy) is the 
potential-flow surface velocity. 

The basic differential Eq. (56) has a singularity at 
¢ = 0. This is evident from the fact that \(dd\/df) 
approaches a finite value and A approaches 0 as ¢ 
approaches 0. Thus, \ behaves as ¢!/? at ¢ = 0, and its 
derivatives are all singular at £ = 0. This presents an 
obstacle for numerical analysis with ¢ = 0 as the start- 
ing point. To overcome this difficulty, the trans- 
formation 


A= (58) 


is introduced into Eq. (56). The series solution for A 
near ¢ = 0 is found to be 


A = a’ — 4C%? + (4/15a)(a1 + (32/a*)C%)g5/? — 
(2/9) [(2a1 + 1)C + (80/a*)C4]g* + O(F7/2) (59) 


where a’ is a constant equal to (dA*/df)-~9. Thus, fora 
given value of a; and a given value of C, Eq. (59) can be 
used to calculate A, dA/df, and d*A/df? at a small ¢ by 
assuming a value of a*. This furnishes the initial data 
needed for step-by-step numerical integration starting 
from this point. The value of a* must be adjusted to 
lead to the correct outer boundary conditions for A. 

Numerical results for the case aj = —1 have been 
obtained by using IBM 650 computers. The values of 
C used are: —0.07071, —0.02828, 0, 0.02, 0.02828, 0.05, 
0.07071, and 0.085. For the numerical integration, the 
starting point is chosen to be at ¢ = 0.0001, and the 
Runge-Kutta method is used. Some discussion of the 
numerical analysis has been given in reference 5. It 
may be mentioned that at least three trials (with three 
values of a”) were used to make an interpolation for a’. 
This value of a? was then used for a check run on the 
computers, and readjusted again if necessary. The 
sensitivity of the velocity profile at the outer boundary 
to a change of a’ near its ‘‘correct’’ value is almost 
linear. For example, for C = 0.085, a 0.05 percent 
change of a? results in approximately 0.025 percent 
change of dA/df at A = 1 (+0.0001). 

For C = 0, Howarth’s numerical results for the 
Blasius profile were obtained accurate to the third 
decimal place, when a? was chosen to be 0.9392 and the 
integration step was taken to be 0.00012. However, in 
order to cut down the machine time, larger integration 
steps were generally used for C # 0. Thus, according 
to Fig. 6, the wall skin-friction coefficient for C = 0 
differs from the Blasius value by about 1.2 percent. An 
error of approximately the same amount is probably 
expected for all the wall skin-friction coefficients shown 
in Fig. 6. Because of the complexity of the problem, 
error estimate is difficult and has not been attempted. 

The velocity profiles v; are plotted in Fig. 3. Also 
shown are the potential-flow velocity profiles computed 
by Murphy’s formula 


u = U(x)/[1 + K(x)y] (60) 


where x and y are the conventional boundary-layer 
coordinates, and A(x) is the surface curvature. The 
difference between v; and u is of the order R~' for R> 1. 

These profiles are similar to those given by Murphy,’ 
but there appear to be some basic differences. Fig. 3 
indicates that the boundary-layer thickness increases 
with C, while Murphy’s results show the opposite 
trend. A close examination of the profiles near the 
wall reveals a rather peculiar behavior, namely, there 
are cross-overs of the profiles as can be seen from Fig. 4. 
In addition, the velocity profiles for C> 0 clearly have 
inflection points away from the wall. Fig. 5 is a replot 
of Fig. 4 by using the abscissa of Fig. 3. No details of 
the profiles near the wall, however, are available from 
reference 2. 
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The wall skin-friction coefficients obtained by the 
present analysis and by Murphy’s analysis.are compared 
in Fig. 6. The difference is rather unexpected. Ac- 
cording to the present analysis, a convex surface ex- 
periences a skin friction above the Blasius value, and 
for a concave surface, the skin friction drops below the 
Blasius value as C increases from zero. But it appears 
to reach a minimum at C ~ 0.5, and increases again 
beyond this point. This feature is also noticeable from 
the velocity distribution shown in Fig. 5. Murphy’s 
results, however, indicate that the wall skin friction for 
a convex surface would be lower than the Blasius value, 
while for a concave wall it is higher. 


Discussion 


The disagreement between Murphy’s results and 
those obtained here is evidently due to the use of 
different boundary-layer approximations and the subse- 
quent different numerical methods of obtaining the 
solutions. It was anticipated at the beginning of this 
work that the wall skin-friction distributions obtained 
from these two different boundary-layer approximations 
might be the same as in the flat-plate case when results 
obtained by Blasius and by Carrier and Lin® are com- 
pared. However, this anticipation turned out to be 
false. 

As far as the boundary-layer equations are con- 
cerned, it is rather difficult to decide in favor of one 
approximation rather than the other. The same 
approximation technique is used, although carried out 
in different coordinate systems. The difference is, of 
course, that in the present analysis no approximation is 
applied to the curvature terms in the equations of mo- 
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Fic. 5. Boundary layer velocity profiles near the wall. 
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tion, Eqs. (5) and (6). Thus, the artificial singularity 
in the differential equation mentioned in the Introduc- 
tion does not appear. The order of magnitude of the 
surface-curvature effects is brought out, however, by 
the boundary condition, Eq. (24) (see the remarks by 
S. Goldstein in reference 10). 

One advantage of the method of approach adopted 
here appears to lie in the boundary conditions. In 
particular, the condition of irrotationality at the ‘‘edge 
of the boundary layer”’ is exact and consistent with the 
vorticity term in the boundary-layer equations. More- 
over, the mathematical problem formulated in the 
generalized coordinates with \ as the dependent variable 
becomes of the same type as the flat-plate problem. 
Thus, it is feasible to obtain satisfactory numerical 
results. However, numerical integration of Eq. (56) is 
by no means an easy task. Additional numerical 
analysis will be done by using the generalized Falkner- 
Skan equation, Eq. (A-3), given in the Appendix. 

Although Murphy’s theory appears reasonable, some 
steps in his method of numerical analysis may be open 
to question. For example, consider the equation 
[Eq. (45) of reference 2] obtained by Murphy for the 
determination of the skin-friction parameter C) (C, = 
1.328 for a flat plate) as a function of the surface- 
curvature parameter A (= —C/2"/*) 

C2?/°[1.65562 — 0.1444 C23? — 1.002An] > 
[2/(1 + 2An)], n>3 (61) 


where 7 is essentially the Blasius variable. Evidently, 
for a given A, the value of C, satisfying Eq. (61) depends 
on the choice of 7. Murphy’s results, A ~ —0.0235 and 
C, ~ 1.52, are obtained when 7 is chosen to be 3. For 
C, = 1.52, if 7 is chosen to be 5, however, A is found to 
be --0.016. Moreover, it is noted that for Cy = 1.52, 
Eq. (61) can also be satisfied by A = +0.224 (approxi- 
mately) when 7 is chosen to be 3. A positive A implies 
that a convex surface would experience a skin friction 
higher than the Blasius value in agreement with the 
results obtained in this paper. 

In the initial phase of the present numerical analysis, 
the values of a? used were obtained from Murphy’s 
results. But it was soon found that these values of a? 
cannot lead to the correct outer boundary conditions 
for the velocity profile. Only after some drastic adjust- 
ments was the wall skin-friction distribution finally 
obtained as a function of the surface curvature as 
shown in Fig. 6. 

The result that the wall skin friction on a convex 
surface is below the Blasius value (for the same 
Reynolds number based on the potential-flow velocity 
at the surface) might be explained as follows. If the 
flow is a potential one, there exists at the (convex) 
surface a positive normal pressure gradient. But, in a 
boundary-layer flow, according to Eq. (6) and the 
boundary-layer approximation, 


Op/On + A*h,(Oh;/On) = O (62) 


It follows that the pressure gradient vanishes at the 
surface. Thus, the boundary layer over a convex 
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Fic. 6. Wall skin friction. 


surface is subjected to a higher pressure as the radius of 
the surface curvature is reduced. This is also apparent 
from the fact that, at the ‘‘edge of the boundary layer,” 
the ‘‘free stream”’ velocity is below that of the potential 
flow at the surface. Consequently, it appears reason- 
able that the boundary-layer thickness should be 
reduced and the wall skin friction increased. On the 
other hand, for a concave surface, the opposite is true. 

Physically it is hard to explain why the skin-friction 
distribution for a concave wall should exhibit the 
curious behavior shown in Fig. 6, namely, the existence 
of a minimum skin friction at C ~ 0.05. No satis- 
factory mathematical proof has been found yet. 
Additional analytical and numerical work is needed 
before an answer to this question can be given. 


Concluding Remarks 


The main purpose of this paper is to present the 
boundary-layer theory proposed by the first author*® to 
account for some surface-curvature effects. The 
numerical example used was chosen for its relative sim- 
plicity, and no undue significance should be placed on 
this flow problem. In fact, some reservations on the 
applicability of boundary-layer theory to this problem 
may be raised. 

As mentioned earlier, the method adopted in this 
paper constitutes an “inverse” approach to the bound- 
ary-layer problem. In this approach, the outer flow is 
not specified beforehand, except that the boundary- 
layer flow is required to approach a potential flow at 
free stream. Thus, the problem of correcting the 
outer flow for the displacement effects, which are 
known to be of the same order of magnitude as the 
surface curvature effects, does not arise. However, if 
the flow problem is of the ‘‘direct’’ type—i.e., if the 
problem is to solve the boundary-layer flow over a 
given body with given “‘free-stream” conditions—the 
outer flow should be corrected with respect to the dis- 
placement effects and also all other possible effects of 
the same order of magnitude—e. g., the free-stream 
vorticity effect—in order to have a consistent analysis in 


order of magnitude. A general analysis of this “direct”’ 


problem appears to be a very difficult task and is not 
yet available.* Two main features of the present 
study by an “inverse”? approach are: (1) it reduces 
the complexity of the problem by isolating the surface- 
curvature effects, and (2) it represents a nonlinear 
analysis of the surface curvature effects. Thus, the 
present analysis may be useful for future work toward a 
more general boundary-layer theory. 

Extension of the present theory to flows with free- 
stream vorticity, to three-dimensional flows, etc., is of 
interest and is being carried out. It is sufficient to 
mention here that for a body of revolution the corre- 


sponding boundary-layer equation of motion is 


P 8 (4)] P ( a 
hava* + Va = 
Oa og he op Oa 


] oO jhig O [ } O (**) | (G2 
. lal 4 2 
R 08 tha 08 656 \ha) If 


where a and £ are the coordinates in the meridian plane 
and §-lines coincide with the streamlines. Let r(a, 8) 
be the metric measure along the circumferential direc- 
tion, and y, be the Stokes stream function defined by 


Va rhg(Ow;/OB) (64) 


= 0 (65) 
If the transverse curvature is ignored—i.e., r(a, 8) = 
r(a)—the following transformation (one form of the 


Mangler transformation!!) 
dé = r—*da,_ dn = r—dB (66) 


reduces Eq. (63) to its two-dimensional counterpart. 
In general, it is desirable to include the transverse- 
curvature effects, and for this case a modified Mangler 


transformation may be used. 
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Appendix 


A Generalized Falkner-Skan Equation 


It is of interest to see that the differential equation 
(39), with H eliminated by using Eq. (40), can be 
transformed to a convenient form which reduces to the 
Falkner-Skan equation® for zero surface curvature. 
Eq. (52) suggests the use of the variables f and x defined 


by 
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i dc 
A = df dx, x = f os (A- 1) 
; oA 


By choosing f = ¢, Eq. (39) becomes 


(e* xf") " (A-2) 


2Cxpr yp 7 2Cx¢ 
aye ry” — fle » ae ‘a 


where the primes denote differentiations with respect 
to x. Eq. (A-2) can be written in the following form, 
by using Eq. (53): 


f+ ff!’ + (1 — 28) ff" + 2C(2f'" + ff" + 
2Cf") = 0 (A-3) 


This differential equation is obviously a generalized 
form of the Falkner-Skan equation. The boundary 
conditions are f = f’ = Oat x = O, and f’ 
as x > ©. This equation, only recently obtained by 
the first author of this paper, was not used, however, in 
obtaining the numerical results presented in this paper. 
Eqs. (39) and (40) were used as the basis of the numer- 


=> 1, f” +0, 


ical analysis.® 





7 7 
Supersonic Shea r Flow (Continued from page 863) 
S(O) = yp: Cri (37) 1= 3, mw = 0.0363, My = 1.35, M,(3) = 1.15 
The solution for \> wis The calculation is straightforward, and only a few 
a eigenvalues were needed to assure convergence of Eq. 
— ee =: . = gt 29 ; : 
n=Ae* sin [(B Vr u* | (38) (11). In Fig. 3 pressure curves are shown. 


With the help of the boundary conditions, Eq. (31), the 
eigenvalues and the constant B are determined as 


A = [w? + (an/r)?)}?, 2 = 1,2,3... (39) 


9 9 


B = (d? — w?)—/? tan—![u-"(A? — gw?)!/?] (40) 
The constant A and the spectrum function S(\) are 
given by substitution of Eq. (38) in Eqs. (34) and (35), 
and by use of Eqs. (39) and (40); we get 


A? = 21 Myer-(My? — 1)-? (41) 
SQ) = 2yprAual—r—nd-3[(— 1) "e"” — 1] (42) 


In the formulas above, 7 is given by Eq. (33), and is 
obtained by substitution of the value of M, at y = 1. 
Now suppose further that the airfoil section is 


parabolic—i.e., 
h(x) = 4hex(1 — x) 

where hy is the maximum thickness of the airfoil. Eq. 
(22) then gives for the function F(x, 0; A): 
— F(x, 0; ) = co cos Ax + (ao/A) sin Ax + co(xA/2)? X 

[1 — (x?d?2/2-41)5 + (x*A4/2-6!)(11/2) — ... ] 
The pressure on the airfoil is finally obtained with Eq. 
(11). Inthe numerical calculation, the following values 


were used 


(5) Discussion 


The pressure curves obtained for the simple flow 
model show a few interesting features. The effect on 
the pressure coefficient in the chord direction of the 
nonuniform flow is a small disturbance of the Ackeret 
two-dimensional airfoil theory. This can be seen from 
Fig. 3, where the actual c, curve at the wall y = 0 is 
compared with its two-dimensional counterpart com- 
puted for the Mach number Mp. The change of c, in 
the span direction is mainly due to the definition of c, 


as a local property. 
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Alignment of Inertial Guidance Systems by 
Gyrocompassing—Linear Theory’ 


ROBERT H. CANNON, JR." 
Stanford Unwersity 


Summary 


The performance of an inertial guidance system can be only 
as good as the accuracy to which the system is initially aligned. 
Optical alignment can be performed with high precision, but 
involves external equipment which may be too elaborate for 
many operational applications. The alternative procedure of 
gyrocompass alignment uses the system’s own instruments, and 
can be accomplished without any external equipment. When 
sufficient alignment time is available navigation errors due to 
gyrocompass alignment are no greater than other errors inherent 
in the system’s instruments 

In certain tactical operations it is desirable to sacrifice some 
accuracy in the interest of more rapid alignment. But, on the 
other hand, if alignment must be performed on a moving base, 
base motions may be filtered only by lengthening alignment time 

In this paper some simple linear gyrocompassing schemes are 
described, and the basic dynamic relations are established for the 
dependency of alignment accuracy upon instrument charac- 
teristics, speed of alignment, and base motions during the align- 
ment period. Upper bounds are established for alignment per- 
formance achievable under typical circumstances. 

It is planned that a subsequent paper will extend the discus- 
sion to describe nonlinear effects and sampled-data techniques, 
and to present operational results including high-latitude, at- 
sea, gyrocompass alignment of inertial ship navigation systems. 


Introduction 


iow PERFORMANCE accuracy of an inertial guidance 
system can be only as good as the accuracy with 
which the system is initially aligned. f 

For the usual Schuler-tuned system for terrestrial 
guidance, an error ¢, in initial vertical alignment of the 
stable platform produces errors in the system’s report of 
distance traveled which vary sinusoidally with time 
according to the relation: 


Ax = R¢@,, (1 


— COS ap) 


in which Ax is the navigation error, R is the earth’s 
radius, and w is the Schuler frequency. (w) = 27/T, 
where 7 is 84 min.) 

For an error in initial azimuth alignment of the plat- 
form, ¢,9, a navigation error is produced which has the 
form: 


Received August 1, 1960. Revised and received July 13, 1961. 

t Based on work done for Autonetics, a Division of North 
American Aviation, Inc.; RCA Airborne Systems Laboratory; 
and Cornell Aeronautical Lab. To be followed by ‘Sampled 
Data Techniques and Operational Performance,’”’ by C. S. Bridge. 
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Ro.o[sin Qet — (Qe/wo) sin wol] 


in which Q, is the earth’s angular velocity. 

These relations are plotted in Fig. | for a particular 
case. That is, if the initial leveling misalignment were 
1 milliradian, then the navigation error would vary 
with time as shown in Fig. la, with a maximum magni- 
tude of 7 nautical miles. If the initial azimuth mis- 
alignment were | mr, then the resulting navigation error 
would look like Fig. 1b, increasing with time initially at 
the average rate of 0.9 nautical mile per hour. (During 
the first few hours of navigation an azimuth misalign- 
ment has the same effect on navigation accuracy as 
drift rate in the leveling gyro of da, = Qz¢.,.) 

For systems which are not Schuler tuned—e.g., in 
vehicles being fired along a ballistic trajectory or into 
orbit—initial errors in alignment are reflected in the 
aiming accuracy of the rocket boost phase and affect the 
ultimate mission accuracy in a complex manner, de- 
pending upon the mission. 

For test flights of guidance systems it is customary, 
prior to flight, to transfer reference direction optically to 
the system’s stable platform from carefully surveyed 
equipment, firmly fixed to the ground. This method is 
certainly the most precise possible. 
operational technique optical alignment is useful only in 
certain instances—e.g., in firing missiles from stationary 


However, as an 


ground installations. 

In many operational situations it is desired to align 
the inertial system without the use of external equip- 
ment, and to do so on a ship or aircraft which is in 
motion. <A typical case is the carrier-deck launching of 
an airplane with an inertial autonavigator in it. The 
airplane’s inertial platform must be aligned while the 
airplane is at any arbitrary location on the pitching, 
heaving deck, and while the carrier is maneuvering. 
(Even for land-based launch, the flexibility gained by 
avoiding external equipment may be of great opera- 
tional importance. Moreover, wind gusts on the field 
may produce serious base motions. ) 

In such a case the problem of inertial system align- 
ment may be solved by employing a ‘“‘gyrocompassing”’ 
This consists, basically, of using the inertial 
platform’s own accelerometers to slave it to the local 
gravity vector, and then letting one of its gyros seek 
true north by sensing the rotation of the earth just as a 
gyrocompass does. The computations involved are 
simple, and are readily incorporated into the inertial 
system’s computer so that no external equipment is re- 


technique. 
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quired. The fundamental difficulties encountered are 
(1) the basic limit imposed on accuracy by gyro drift 
and inaccurate velocity information, (2) the sensitivity 
to base motions of a system designed for very fast 
alignment, and (3) the ineffectiveness of the method at 
very high latitudes. 

Other situations in which the gyrocompass technique 
is used include preflight alignment of long-range, air- 
launched missiles, and the alignment of ship-borne 
inertial autonavigators when land and celestial fixes are 
unavailable. Gyrocompassing was used extensively, 
for example, to align the inertial systems during the 
subpolar voyages of the submarines Nautilus and Skate. 

More generally, gyrocompassing may be employed 
whenever a platform or vehicle is to be controlled to a 
reference which is rotating about a known axis, as in a 
local-level satellite. 

The accuracy with which an inertial navigator can be 
aligned by gyrocompassing is determined, basically, by 
the accuracy of the system’s own inertial instruments— 
the accelerometers which are used to locate the gravity 
vector, and the gyro which is used to sense the earth’s 
rotation—and by the accuracy with which base velocity 
is known. Thus, given enough time, and accurate 
velocity information, gyrocompass alignment accuracy 
can be made comparable with the inherent navigation 
accuracy of the system. 

In some operations ample time is available for align- 
ment: on shipboard it is reasonable to allot several 
hours at a time; missiles and aircraft at ‘‘power-on 
standby alert’? may keep their systems in continuous 
alignment. 

But in certain other operations speed of alignment is 
all-important—as, for example, when aircraft must be 
launched from a power-off status in a matter of minutes. 
In such cases a trade-off must be made between align- 
ment speed and system accuracy. In the case of align- 
ment on a moving base, the motion characteristics of the 
base are a key parameter in this trade-off. 

In this paper some simple, linear gyrocompassing 
systems will be described and the basic dynamic rela- 
tions between accuracy and instrument characteristics, 
speed of alignment, and base motion will be derived. 
These relations represent fundamental limitations. 

Before discussing any alignment system the general 
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geometric and dynamic properties of an inertial plat- 
form will be derived. 

Then, in order to develop the major alignment rela- 
tions in the simplest mathematical terms, a linear, 
second-order control system will be described for 
effecting gyrocompass alignment, and its capabilities 
will be analyzed. 

After that a more practical, third-order control 
system will be investigated. 

Finally, the importance of nonlinearities and of 
digital components will be described. 

In a subsequent paper by C. S. Bridge the discussion 
will be extended to describe more fully the effects of non- 
linearities and the gains attending the application of 
sampled data techniques, and to present operational 
results including high-latitude, at-sea gyrocompass 
alignment of inertial guidance systems on ships and 
submarines. 


Platform Dynamic Relations 


Alignment Geometry 

Fig. 2 shows a typical stable platform, controlled by 
single-axis gyros—gyros whose spin axes have only one 
degree of freedom with respect to the platform. (Con- 
trol by two-axis gyros is also common and results in the 
same alignment dynamics. Single-axis gyros are chosen 
for illustration here because the geometery is somewhat 
easier to follow.) 

Each gyro in Fig. 2 controls one axis of the platform 
via an appropriate system of gimbal servos. For 
example, the z gyro in Fig. 2 senses’ any azimuth mo- 
tion of the platform and informs the z-axis controller, an 
electrical system which in turn applies a corrective 
torque via the z gimbal motor. The x and y gyros 
control outer gimbals which have been omitted from 
the picture. 


* By precessing about its output axis. 

™ Cannon, R. H., and Chandler, D. P., Stable Platforms for 
High Performance Aircraft, Aero. Eng. Review, Vol. 16, No. 12, 
Dec. 1957. . 


Z-AXIS CONTROLLER 











Fic. 2. Stable platform, using single-axis gyros. 
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The time of response of each platform servo—involv- 
ing time to accelerate the inertia of the platform, ete.— 
is commonly shorter than 0.01 sec. Since alignment 
time will be measured in minutes we can safely neglect 
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Fic. 3. Coordinates for gyrocompassing study. 


servo response time as a factor in our study of alignment 
dynamics and assume that the platform is always exactly 
aligned with the input axes of its gyros. These are the 
x-, y-, and z-axes in Fig. 2—that is, the x gyro senses 
motion about the x-axis of the platform, etc. By the 
above assumption, x, y, z constitute a platform-fixed, 
orthogonal coordinate system. 

The control system is supposed to effect alignment to 
the local vertical and to north. It can be shown 
(Appendix A) that alignment about the azimuth and 
east-west axes of the platform is tightly coupled (due to 
earth rotation), but that alignment about the north- 
south axis is almost independent—that is, of the desired 
alignment is z-axis vertical, x-axis north, in Fig. 2, then 
alignment about x is almost independent of alignment 
about y and g. Furthermore, alignment about x 
involves seeking the vertical only (no gyrocompassing) 
and can be accomplished very rapidly. It is therefore 
reasonable to assume, for the purposes of the present 
study, that x-axis alignment has already been accom- 
plished and is held firmly during subsequent y- and z-axis 
alignment. This assumption will greatly simplify the 
study without altering its conclusions. 

Fig. 3 is drawn to define carefully the alignment 
geometry. First, an earth-fixed coordinate sphere 
(“unit sphere’’) is imagined to be drawn with its center 
at the center of the platform (and also of the y gyro). 
On the sphere are painted the earth-fixed coordinates 
NWV: N pointing north, W west, and V vertically up- 
ward (Fig. 3a). Next, points x, y, z are marked where 
the axes of the platform pierce the sphere. When 
these coincide with NWV the platform is aligned. 
When they do not, azimuth misalignment is defined 
simply as the angle, $., between vertical planes xV and 


NV; leveling misalignment ts the ‘‘dip’” angle, dy, be- 
ween the x-axts of the platform and the horizontal plane, 
NW. (¢,; is zero.) 


Gyro Dynamics 

The key dynamic relations in Fig. 2 are the torque 
balances on the output axis of each gyro. These con- 
sist of a gyroscopic torque, an uncertainty torque, and a 
control torque from the alignment system. 

The gyroscopic torque is /J XQ, in which H is the 
spin angular momentum of the gyro, and Q is the total 
angular velocity of the platform, including earth rota- 
tion, Qg, and ship motion, Vy/R and —Vy/R. In 
Fig. 3a, Qg has been resolved into its vertical and 
horizontal components through the latitude angle, X. 

Uncertainty torques, U, include such unwanted 
effects as unbalance, friction, stray magnetic fields, 
temperature effects, etc. These are, of course, made 
as small as possible, and the figure of merit for a gyro- 
scope is the drift rate it exhibits due to uncertainty 
torques: 


2, = U/H 


seu 
Control torques, C, are applied via electromagnetic 
torquers on each gyro output axis. During normal 
navigation these torquers are used to compensate for 
earth rate and for steady, known bias torques. During 
alignment additional electric current is fed the torquer 
from the alignment system to control platform motions. 
Referring to Figs. 2 and 3b, and ignoring Vw for the 
present, the angular velocity of the platform about its 
y-axis is (for small angles ¢, and ¢,): 


dy — o:Q¢ cos X + (Vx/R) 
and the gyroscopic torque about the y-gyro output axis 
is 
H(¢, — $2 cos X + (Vy/R)] 
so that the torque balance about the y-gyro output axis 
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4(b) MATHEMATICAL REPRESENTATION 


Fic. 4. Simple gyrocompassing system. 
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H¢, — H¢Mez cos \ + H (Vy/R) — U, — C, = 0 (1) 


in which U, means “uncertainty torque about output 
axis of y gyro” and is considered positive when it pro- 
duces positive drift rate, ¢d,. In operation the best 
available independent velocity measurement is used to 
bias the y gyro to compensate for Vy/R. The error in 
this corrective bias, Vy,/R, may also be considered 
partof Uy. This is discussed further below. 
Similarly, for the z gyro: 


H(¢: + Qe sind + ¢,Qz¢ cos A) — U, — C,’ = 0 


or, with fixed earth-rate compensation included in the 
control torque* 


(C,! = C, + AHQg sind): 
Hd, + Hd Qe cos \ — U, — C, = 0 


These expressions are more conveniently written in 
terms of equivalent gyro rates, 2, = U/H,Q, = C/H: 


dy = 282, + Quy = Qey 


: ina — $2, + Qus + 2. (la) 


in which Q, is ‘‘local earth rate’: Q, = Qg cos Xd (Fig. 
3b). (Again, errors in compensation for Vy/R and for 
earth rate on the z gyro are included in Q,, and Q,,, 
respectively.) 


‘‘North Steaming”’ Error 


In the development above, only error due to north- 
ward vehicle velocity, Vy, is considered. In general, of 
course, the vehicle may be traveling with velocity Vw 
west, as well as I’y north, so that the platform angular 
velocity will have both additional components — Vw/R 
and Vy/R = about the N and W axes, respectively, 
in Fig. 3 (because the platform is slaved to the local 
vertical). 

The total effect of component Vyj/R is simply to re- 
duce the magnitude of multiplier 2, in Eq. (1). Ona 
ship the effect is not significant (Vw/R)/Q, being of 
order 0.1 at most (except at very high latitudes). But 


* Again, any error in this compensation—e.g., due to torquer 
anomalies—becomes part of the uncertainty torque, U,,. 
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Fic. 5. Dynamic characteristics of system of Fig. 4. s(s + Ky) 
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Fic. 6. Response of system of 
Fig. 3 to vehicle motions. 








in a missile, of course, gyrocompassing may be en- 
hanced in an easterly flight, or greatly weakened in a 
westerly flight. (In particular, a missile which flies at 
minus local earth rate has no inertial angular velocity 
and gyrocompassing is impossible.) In such cases the 
same analysis may be used merely by replacing 2, by 
Q, — (Vw/R)| everywhere. 

The effect of component Vy/X, on the other hand, is 
to change the apparent direction of the 2, vector, which 
is always critical if Vy is not zero. This effect must be 
compensated for, as it must in conventional gyro- 
‘north 


compasses where the compensation is known as 
steaming error correction.’’ The best available inde- 
pendent measurement of north velocity is used—e.g., a 
pitometer or Doppler radar, and a bias torque, corre- 
sponding toa platform rate —\ = — Vx/R, is applied to 
the y gyro (Fig. 2). Any error in this bias signal—i.e., 
in velocity information—is part of Q,,, and will have 
exactly the same effect as a y-gyro drift rate. That is, 
an error of 1 knot in north-velocity information has the 
same effect as a y-gyro drift of '/s0 degree per hour. 


Second-Order Control System 


A simple scheme for accomplishing both leveling 
about the east-west axis and azimuth alignment is 
shown in Fig. 4a. The scheme uses the system’s own 
instruments: leveling gyro, azimuth gyro, and ac- 
celerometer. The torque balances shown (heavy 
lines) on the two gyros represent Eqs. (1)." In the 
absence of disturbing motions, bias, and uncertainty 
torques, the accelerometer output, a, indicates directly 
the platform dip angle, ¢,. (Only angles less than 6° 
are considered. @, isin radians, and a in fractions of a g.) 

To accomplish leveling, the y gyro is torqued with a 
signal proportional to a (K,’ loop in Fig. 4). In the 
presence of an azimuth angle ¢, the leveling gyro 
senses some earth rate so that the level loop will not, in 
general, drive the platform to ¢, = 0. Instead a small 
dip angle will remain to furnish a control torque to the 
gyro just equal to the earth rate torque: 


*The torque —H¢, in the second of Eqs. (1) is not mech- 
anized in Fig. 4. It is omitted because the control torque, 
Tz = —Kz'¢y (Fig. 4), will always be very much larger for useful 
values of K,.’. 
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To accomplish azimuth alignment a torque propor- 
tional to a is applied also to the azimuth gyro as shown 
(K,’ loop in Fig. 4a). Then in steady state d, = 0, ¢, = 
0a and, therefore, ¢, must equal zero. 

In practice the accelerometer may be subject to a 
bias, ad. The accelerometer will also record actual 
lateral accelerations, a,;, of the vehicle. These will 
contribute to the dynamic and static response of the 
alignment system. 

We shall consider, in turn, (1) the dynamic charac- 
teristics of the system of Fig. 4a, (2) the static errors 
which will result from bias and uncertainty torques, (3) 
the system’s response to vehicle motions during the 
alignment operation, and (4) the speed with which the 
system will reach a desired alignment accuracy. Sys- 
tem characteristics will be considered before static 
errors because the magnitudes of the static errors de- 
pend upon system gains K,’ and K,’. Response to 
vehicle motions will be considered before transient re- 
sponse because the system time constant will depend 
upon the amount of filtering to be done. 


Dynamic Characteristics 

To expedite the dynamic studies Fig. 4b has been 
drawn to depict a Laplace transformation of the equa- 
tions represented in Fig. 4a.* Some multiplying and 
dividing by H has been done to simplify the picture. 
Note that K, = K,’/H and K, = K,’/H and that un- 
certainty torques U have been replaced by corre- 
sponding gyro drift rates Q, = U/H. The Laplace 
transform of these steady drifts is Q,/s. Capital letters 
are used to represent the Laplace transforms of corre- 
sponding lower-case letter variables. Initial conditions 
have also been indicated (as impulses) in Fig. 4b. A 
constant accelerometer bias is represented by 4a)/s. 

The characteristic equation of the system of Fig. 4b is 


* Alternatively the ‘‘transfer functions’’ in Fig. 4b. could be 


obtained from Eqs. (1) by assuming solutions of the form ¢y = 
$, e*', etc. 
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s(s+ K,) + KQ, = 0 (3) 


The relation between gains A, and A, and the charac- 
teristics of the system is indicated in Fig. 5, in the ‘‘s 
plane,”’ which plots vertically the imaginary parts (fre- 
quencies) of the roots of Eq. (3), and horizontally the 
real parts (decay rates). It is clear that the charac- 
teristic roots of the system can be given any value 
desired by appropriate choice of K, and A,. 

It will be found that the best transient response will 
be obtained by using characteristics somewhere between 
those labelled s; and s:in Fig. 5. Any such roots will be 
characterized by a damping time constant r where, 
from Fig. 5, 7 is determined by the choice of K,: 


K, = 2/r 
and, for roots s; = —1/r, —1/r, K; is given by: 
KQ. = 1/7?) 
ar : (4) 
So that: K,/K, = 2rQ, 
(Roots ss = —1/r + 7 (1/7) would be obtained by in- 


creasing K, to the value 2/7°2,. It will be found that 
the difference in performance between roots s; and Se 
is not of major importance in this study. Roots s, will 
be used henceforth for ease of mathematical manipula- 
tion.) 

If we take a nominal value of earth rate, Q, = 5 X 
10—* rad/sec (at 45° latitude), relations (4) become: 


K, = 2/r 
K, = (2/r)(104/r) ? (4a) 
K,/K, = 1/108 


Also, the open loop relationship (2) between ¢, and ¢, 
becomes: 


¢. = (K,/Q.) dy = (4 X 10~4/7)¢, (2a) 
Static Errors 

‘ . . gB.. ° 

Static errors can be obtained most directly from the 
torque balances shown in Fig. 4a. In the presence of 
gyro uncertainty torques and accelerometer bias, these 


are: 
(a) HQ.¢. + U, + K,’a = 0 
(b) + K,a+ U,=0 





8b 


Fic. 8. Elements of transient response. 
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Or, combining (a) and (b) and dividing by HQ,: 
pz — (Qyy/’ 2.) + (Ky, ‘K = (Quz 2.) / (5) 
» (5 
and dividing (c) by K,’: @ = @ + (uz K.)\ 


For rapid alignment Eqs. (5) can be simplified because 
the last term in each turns out to be very small. For 
example, consider a system time constant rt = 100 sec. 
Then, from Eqs. (4a), Kz = 2 and Ky/Kz = 0.01. 
Then in the first of Eqs. (5) the z-gyro drift could be 10 
times as great as y-gyro drift and still be relatively unim- 
portant. In the second of Eqs. (5) an azimuth gyro 
drift as great as 1° per hour would be comparable to an 
accelerometer bias of only 1 part in 200,000. Thus, for 
a fast-alignment system Eqs. (5) can be approximated 
by: 


¢ = —- (Quy/Q.) = Qya + Vue /R)/S Q, 


in which Q, is drift rate of the y gyro, and Vy, 


(5a 
is error in knowledge of north velocity; and by{ ~* ) 


gy = % 

For example, to obtain 1 mr alignment accuracy in 
azimuth at 45° latitude the leveling gyro drift rate 
must be less than 0.01°/hour, and north velocity must 
be known to 1 ft/sec. To obtain 1 mr leveling accuracy 
ad must be less than 10~*g. 

Since Qe = Qe cos X — Vyw/R, the difficulty of 
accurate gyrocompassing at very high latitudes is 
apparent. 


Response to Vehicle Motions 

The response of the system to lateral accelerations of 
the vehicle in which it is mounted can be obtained in 
transform form directly from Fig. 4b: 











b,(s) =< Ke + Km _2rs +1 
j [s(s + K,) + KQ)°% (s+)? 
®,(s) = ~ ————~ Ay 

[s(s + Ky) + Ka] 4 ” s+ 1)? 


The expressions to the right result from using values of 
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the K’s given by Eqs. (4). If the disturbing motions 
are considered to be sinusoidal, of the form a; cos w/, 
the peak magnitudes are given by: 


¢ | 2jwyT + 1 | 
'y|mar 
) ] ‘ 
Gay + 4 (6) 
| jos Q, 
gz maz — - , ay 
(Jrw; + 1)° 


These relations are plotted asymptotically in Fig. 6. 
It is noted that the azimuth response to lateral disturb- 
ing motions is very much greater than the leveling re- 
sponse for typical values of 7, the time constant of the 
alignment system. 

For design purposes it is convenient to redraw Fig. 6 
as a plot of ‘‘allowable motions.’’ Consider, for 
example, the nominal case that azimuth swings due to 
vehicle motions are to be maintained below 1 mr. 
Then Fig. 6 can be replotted as in Fig. 7 where the solid 
lines indicate the maximum allowable lateral accelera- 
tions of the vehicle during alignment with a second- 
order system. 

Some typical amplitudes of disturbance motion have 
been plotted as points on Fig. 7 for comparison. The 
ship motions, plotted as single sine waves, correspond to 
the natural rolling of a ship. The frequency chosen for 
the airplane motions corresponds to a rigid body 
natural frequency of the airplane on its tires in the 
presence of gusty wind. (A more sophisticated ap- 
proach would involve statistical analyses of gust spec- 
tra.) 

It is seen that for this system, alignment on heavy air- 
craft can be accomplished with fast time constants (r > 
100), but that for alignment on light aircraft, or on ship- 
board, longer time constants (r & 1,000) must be used 
or else the vehicle must be sheltered from wind and 
waves. 

A more complex situation, for both ships and air- 
craft, involves vehicle maneuvers under way. For 
alignment at sea or in flight it is necessary to bias the 
accelerometer for vehicle motions, using as a reference 
some other measurement of motion—e.g., a pitometer or 
radar system. The disturbances to be plotted in Fig. 7 
would then consist of the difference between actual and 
estimated motions. In this case the magnitude of 
allowable motion will depend on the accuracy with 
which vehicle motions are known during alignment. 
In particular, certain very slow maneuvers may catch 
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Fic. 10. Third-order gyro compassing system. 
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Fic. 11. Dynamic characteristics of third order systems. 


the curve in Fig. 7 at its lowest point where—for 7 = 
100, for example—allowable accelerations are in the 
micro-g class. 

From Fig. 7 it is clear that the simple system of Fig. 
4 is severely hampered in that quite long time constants 
may be required to perform necessary filtering of 
vehicle motion. As an alternative to using long time 
constants, the order of the system can be increased to 
provide steeper filter lines in Fig. 7. Design of a third- 
order alignment system will be discussed in a later sec- 
tion. 

First, we shall complete our study of the simple sys- 
tem of Fig. 4 by computing its transient response, and 
then its total alignment time. In these calculations 
typical r’s are to be chosen on the basis of Fig. 7. 
Transient Response 

The transient response to initial values of ¢, and @¢, 
can be written in Laplace form directly from Fig. 4b: 
7 TS$,(0) + 7Q.6-(0) 

(rs + 1)? ' 


Sdy(0) + 26-(0) 
[s(s + Ky) + K.Q] — 
_ (s+ K,)$.(0) — K.4,(0) _ 
[sis + Ky) + KG) 
(rs + 2) ¢-(0) as (1, 72), (0) 
(rs 1)? . 


®,(s) = 


®.(s) 


Again, the expressions to the right result from using 
gains (4). 

The corresponding transient motions are: 
$y = $y(0) [1 — (t/7)] e7”" + @.(0) (72) (t/r)e"" (7a) 
, = $,(0) [1 + (t/r)] e~ 7 —,(0) (1/72) (t/r)e~””_ (7b) 


The coefficient 72,, which multiplies ¢,(0) in Eq. 
(7a), is usually a very small quantity. Remembering 
that at 45° latitude Q, is about 1/(2 X 10‘), this 
term can be neglected for reasonable values of r—say, 
7 < 2,000 sec. 

On the other hand, the term 72, appears in the 
denominator in Eq. (7b) so that the initial level angle, 
¢,(0), has a preponderant effect on the azimuth tran- 
sient. For this reason it will be found expedient to per- 


form a preliminary leveling operation before closing the 
azimuth (K,) loop in Fig. 4. Preliminary leveling will 
be discussed below. 

The important transients from Eqs. (7) are plotted 
in Fig. 8. In particular, the azimuth transient will 
consist of some combination of curves B and B;. These 
curves are shown on the usual linear scale in Fig. 8b. 
But since we shall be interested in very small final 
magnitudes which may not be reached for eight or ten 
time constants, the curves are also plotted with ampli- 
tude to a logarithmic scale in Fig. Sa. 

(The effect of using equal real roots from Fig. 5, 
rather than a complex pair, can be judged from Fig. 8. 
A system having 0.7 critical damping will oscillate be- 
tween limits whose magnitude is V2 as large as curve 
A, Fig. 8. Thus, for a typical transient range of, say, 1 
to 0.01, the critically damped system would take 6'/2 
time constants (curve B;) while the 0.7 damped system 
would take only 5. But it must be remembered that 
the choice of 7 will depend upon the filtering required of 
the system and that to obtain equally effective filtering, 
the system having ¢ = 0.7 will have to have a time 
constant which is \2 times as large as the system 
having ¢ = 1. Thus, the real transient time would 
actually be slightly /onger for ¢ = 0.7.) 

Once the value of 7 has been decided upon, the time 
for closed loop alignment to a specified accuracy can 
be read from curve B, in Fig. 8. 

Some additional time will be required for preliminary 
leveling. This will now be computed. 


Preliminary Leveling 

The object of preliminary leveling is to bring the 
second term of Eq. (7b) down to the same level as the 
first term. The corresponding value of ¢,(0) depends 
of course, on the value of r to be used during subsequent 
closed-azimuth-loop operation. Assuming that ¢,(0) 
and ¢,(0) are about the same size—e.g., 100 mr— 
preliminary leveling must reduce ¢, by the factor: 


(dy) pre. = (7Q,) ¢,(0) 


for example, [r/(2 X 10~4)]¢,(0) (8) 


also, (1) pre. © (72)¢-(0) 


where 7 is the value to be used in subsequent closed- 
loop alignment—e.g., for 7 = 1,000 sec, ¢,(0) = 100 
mr, (y] pre. = 5 mr). 

Preliminary leveling will be performed using the 
system of Fig. 4 with Kz = 0; that is, the leveling loop 
is closed but the azimuth angle is left unchanged at its 
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initial value ¢,(0). The gain, A,,, used during pre- 
liminary leveling determines the system time constant: 
7o = 1/K,,. 

7> must be chosen so that lateral accelerations, a,;, of 
the vehicle do not produce ¢,’s larger than those of Eqs. 
(8) during preliminary alignment. The preliminary 
leveling system obtained by setting AK, = 0 in Fig. 4b 
acts as a simple first-order filter to a,: 


.. == A s/(ros + 1) 


so that the value of 79 required to filter a given disturb- 
ance of amplitude, a;and frequency, wy, is: 


7 = [ay (dy) pre] 1/wy (9) 


The transient response of the preliminary leveling 


system is: 
by = o,(O)e~”” + (72.)¢-(0) [1 — e”] 


in which 7» will be found [from Eq. (9)] to be quite 
short, so that the second term can be neglected : 


$, = ¢,(0)e~’” (10) 


The transient is given by curve A in Fig. 8. 

(It is noted in passing that in preliminary leveling 
¢, does not go ultimately to zero but goes instead 
[from Eq. (8)] to the value ¢, = +702.¢.(0). How- 
ever, this does not reduce the effectiveness of prelimi- 
nary leveling: on the contrary, it can be shown that if 
¢, has this value rather than zero, the subsequent azi- 
muth alignment transient may actually be speeded up 


somewhat.) 


Computation of Total Alignment Time 


Total time to accomplish alignment under a given set 
of conditions is computed in the following steps: 

(1) The system time constant is read from Fig. 7 
for the alignment accuracy desired and for the given 
magnitude and frequency of vehicle disturbances to be 
expected during alignment. 

(2) The number of time constants required to re- 
duce the misalignment from its expected initial value to 
the desired final accuracy is read from curve B, in Fig. 
8. This is then multiplied by the results of step (1) to 
obtain the time for closed-loop alignment. 

(3) The preliminary leveling accuracy required is 
determined from Eqs. (8). 

(4) The preliminary leveling time constant is com- 
puted from Eq. (9). 

(5) The number of time constants required for pre- 
liminary leveling is read from curve A, Fig. 8. This is 
then multiplied by the results of step (4) to obtain the 
time for preliminary leveling. 

(6) The results of steps (3) and (5) are added to ob- 
tain total alignment time. 

A sample calculation is given in Appendix B. The 
relation between required accuracy and time-to-align 
for a second-order system is given in Fig. 9 for various 
levels of vehicle motion. For reasonable accuracy 
levels—e.g., 1 mr—the alignment time is seen to be ex- 
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tremely long. This is because of the poor filtering 
capability of a second-order system (Fig. 7). 

In the next section it will be found that more efficient 
filtering can be obtained by using a third-order system 
so that, even though more time constants are required 
for alignment, the total alignment time will be found 
to be much shorter. 


Third-Order System 


Mathematically a third-order system can be ob- 
tained by filtering the accelerometer signal as shown in 
Fig. 10 (Fig. 10 corresponds to Fig. 4b). It is shown in 
Appendix C that a filter of this type already exists, in 
mechanical form, for systems which use a velocity 
meter instead of an accelerometer. 

The characteristic equation for the system of Fig. 10 is 


s[s(71 5 +1) + K,] + K.Q, = 0 (11) 


The roots of this characteristic equation are plotted in 
Fig. 11, where it is shown that a set of three equal roots 


1/7 = 3/7 | 
K, =1/r (12) 
and Kk Q. = 1/3" | 


(This assumes 7; can be varied as desired.) 

The locus of roots of the K, loop is plotted dotted in 
Fig. 11 as a function of A,. (If asmaller K, had been 
used roots s; would have been nearer the real axis in Fig. 
11 and the K, locus would have gone through a region of 
unequal real roots. If, on the other hand, A, had been 
larger so that roots s; were farther from the real axis, 
then the A, locus would not have produced any real 
roots but would have merely dipped toward the real axis 
and then moved away from it, to the right.) 

The static errors of the third-order system are the 
same as for the second-order system: that is, they are 
given by Eqs. (5). This is because the gyro torque 
balances in Fig. 10 are the same as those in Fig. 4. 

The response of the third-order system to vehicle 
motions is, of course, expected to be considerably 
smaller than for the second-order system. The re- 
sponse relation for ¢, is obtained in Laplace form from 
Fig. 10: 
®, Ks 5/Q 13 
A, —_ s[s(ms +1) +K,)+K2Q, (rs+1)8 — 
[The second form is obtained by substituting from Eqs. 
(12).] Again, if the vehicle lateral motion is sinusoidal 
of the form a; cos wt then the corresponding amplitude 
of platform azimuth motion is: 

ie /Q 
D i eee (13a) 


ay  \(jroy+ 1)3 


As before, it is convenient to plot this relation in terms 
of allowable vehicle motion for a given desired accuracy 
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of ¢., using t (the time constant of the alignment sys- 
tem) as a parameter. The result is plotted as a set of 
dotted lines in Fig. 7. 

The great difference in the value of 7 required to doa 
given job of filtering is evident. Typically, for the 
estimated motion of a light airplane the second-order sys- 
tem requires a 7 of about 500 sec, while the third-order 
system requires only t = 40 sec. This difference is, of 
course, greater for larger disturbances. 


Transient Response 


The Laplace transform for the azimuth transient 
response of the system of Fig. 10 is: 
#(s) = Slus + D + Krld.(0) — Kidy(O) 
7 s[s(ms + 1) + K,] + K.Q, 


from which [using Eqs. (12)] the transient motion is: 


4-00 [1 +£42(4) Jem 
1 |; (‘)] e7/* (14) 


The elements of these transients are contained in the C 
curves in Fig. 8, from which it is seen that a somewhat 
greater number of time constants is required to reduce 
the transient motions to a given level than were re- 
quired by the second-order system (B curves). How- 
ever, the length of each time constant 7 has just been 
shown to be many times smaller for the third-order sys- 
tem, for typical filtering requirements, so that a large 
net gain is realized by using the third-order system (see 
below). 

The factor 72, is seen to magnify greatly the impor- 
tance of initial tilt angle, ¢,(0), in Eq. (14) in the same 
manner as it did in Eq. (7b). For this reason pre- 
liminary leveling is again required. The system in use 
during preliminary leveling is the system of Fig. 10 


with K, = 0. This system has a response to sinusoidal 
vehicle motion given by: 

®,(s) = — A,(s) 1/(tos + 1)? (15) 
where K,, = 1/27. and 1 = 1/2 


and a transient response given by: 
by = $y(0) [1 + (t/r0) ]e-”” 


which is curve B,; in Fig. 8. Again, 7» can be made 
much shorter than the + to be used after the K, loop is 
closed. From Eq. (14) the preliminary leveling opera- 
tion must bring the value of ¢, down to: 


(dy) pre./Py(O) = 72, (17) 


so that it will have no more importance than ¢,(0). 
For sinusoidal disturbance, ay sinw,t, the preliminary 
leveling time constant must be [from Eq. (15)] : 


To = 1/ey [4s/($y) pre} *” (18) 


The computation of total alignment time for a third- 
order system can now be made using the same six steps 


that were outlined for second-order systems, with the 
following changes: 

(1) Read system time constant from dotted curves 
in Fig. 7. 

(2) Find number of time constants required from 
curve C, in Fig. 8. 

(3) Determine preliminary leveling accuracy re- 
required from Eq. (17). 

(4) Determine 7 from Eq. (18). 

(5) Read number of time constants required for pre- 
liminary leveling from curve By, Fig. 8. 

Values of alignment time are plotted in Fig. 9 for 
comparison with the corresponding numbers for a 
second-order system. It is seen that the third-order 
system results in alignment times which are better, 
typically, by an order of magnitude than those of the 
second-order system. 

The question may be raised whether a system of 
greater than third-order might result in even more im- 
proved speed of alignment. The problem in going to 
systems of higher order is a hardware problem which will 
be discussed in the next section, after some critical non- 
linearities have been mentioned. 


Alignment With Arbitrary Azimuth Angle 


It is sometimes possible to save some alignment time 
by accepting, in the autonavigator, a precisely-known 
azimuth misalignment. In this scheme the gyro- 
compass technique is used, not to align the platform to 
north, but instead to measure what the azimuth mis- 
alignment is. This information is then used by the 
guidance computer to resolve the velocity signals (from 
platform-mounted instruments) into their north-south 
and east-west components. 

Fig. 12 shows an alignment scheme in which the azi- 
muth angle of the platform is left arbitrary (presum- 
ably within a few degrees of north by rough alignment of 
the gimbals). The output quantity for this system is 
¢:,, the indicated value of azimuth angle, which is fed to 
the computer to be used for resolution purposes. It 
will be now shown that if it is linear, the dynamic 
response of such a system is the same as that of Fig. 4b, 
or, if the additional filter is used, the same as Fig. 10. 

The error in this system may be considered to be the 
amount by which indicated azimuth angle differs from 
the actual angle. From Fig. 12 the system response is 


given by 


7 (rs + 1)? (ros + 1) 





_ Bet (r(0)/M) + (8/%) Ay ay 


If this expression is compared with Eq. (13), it is seen 
that the filtering characteristics of the system of Fig. 12 
are the same as those of the system of Fig. 10. The 
extra filter in Fig. 12 is highly desirable for the same 
reason as in Fig. 10. By comparing Eq. (19) with the 
equation preceding Eq. (14), it is concluded that the 
transient response of the system of Fig. 12 is also the 
same as that of Fig. 10 and is given by curve C, Fig. 8. 
Thus, in the absence of limiting in the instruments, 
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there is no speed advantage in the arbitrary azimuth 
system. 

The advantage of this system accrues when limiting 
of the z gyro torquer becomes an important factor. 
Then the system of Fig. 12 may eliminate several 
minutes of slewing time (see next section). 


Nonlinearities 


A number of nonlinear characteristics may exist in an 
inertial system which lead to modification of the results 
obtained above. Three of these—gyro torquer limit- 
ing, digital resolution, and sampled data techniques— 
are discussed briefly below. A much more detailed 
treatment of these effects is given in Part II. 

Gyro Torquer Limiting 

If the torque limit of the gyro torquers is exceeded 
during alignment, the time to align will be increased be- 
cause the steep portions of the transient response 
curve—e.g., Fig. Sb—will be made less steep. 

Curve C, of Fig. 8b is reproduced (solid curve) in 
Fig. 13. We want to see what changes would take place 
in this transient if the gyro torquer is limited to some 
particular maximum value. We do this by drawing a 
slope on Fig. 13 corresponding to the maximum gyro 
slewing rate (at the maximum torquer capacity). 
We then can see what areas of the solid curve exceed 
this maximum slewing rate—e.g., the part of the curve 
between points A and B. With the torquer limited, 
the transient will proceed from ¢,4 to ¢g at only the 
maximum available velocity. This will stretch out the 
curve. A convenient way to show this is to plot back- 
ward from point B (dotted curve in Fig. 13) until ¢, 
isreached. Then the rest of the curve can be plotted as 
before. The additional time is available as a shift in 
the origin, At. Such a shift may involve, typically, two 
or three 7’s and thus lengthen the total alignment time 
accordingly. 

Digital Resolution 

The most severe resolution problem anticipated in- 
volves the digital counter on the output of the velocity 
meter. Since a certain change in indicated velocity 
must take place before the meter clicks it is possible for 
the platform to oscillate back and forth through similar 
limits, continuing its motion in one direction until 
enough velocity has been stored to click the counter by 
one step, which then drives the platform back in the 
opposite direction. 

The amplitude of the resulting limit cycle can be 
calculated only by using a trial and error technique. 


Sampled Data Techniques 

A typical sampled data scheme involves first torquing 
the leveling gyro with an impulse for each bit of accel- 
eration indicated by the velocity meter: this portion of 
the operation is analytically similar to the closed-loop 
part of Fig. 12. 

But in addition, the output of the velocity meter is 
stored in a register for a period of time 7’, at the end of 
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Fic. 13. Increase in alignment time due to gyro torquer limit. 


which the total indicated change in velocity is divided 
by 7 and the average value is used by the computer to 
compute current azimuth angle ¢,;. It can be shown 
that the process of averaging with the register is 
approximately equivalent to continuous filtering with a 
first-order filter having time constant 7. (The proof of 
this involves the fact that for a sinusoidal motion the 
averaged value will be exact within + half a sine wave. 


Since the residual half sine wave is averaged over time 7 


its importance is diminished as 7 is increased just as it 
would be for a first-order filter.) 

Thus the sampled data system is approximately 
equivalent to the system of Fig. 12 for the case of no 
limiting, and the performance can be expected to be the 
same for either a continuous or sampled data system. 
In particular the benefit to be gained by an additional 
filter is as shown by Fig. 9. 

In the case of a system with torquer limits the 
penalty paid for slewing limitations turns out to be less 
with the sampled data system. However, as has been 
mentioned, the digital resolution problem may be 
serious. 

Additional Filtering 

The question might be raised whether even more 
filtering than that shown in Fig. 10 might improve sys- 
tem performance beyond that indicated by the dotted 
lines in Fig. 9. 

In certain digital systems one additional order of 
filtering may be readily obtained. These will be dis- 
cussed in the seque! paper by C. S. Bridge. The at- 
tendant benefits are clear by extrapolation of Figs. 7 and 
9. 

For an analog system additional filtering with long 
time constants required—e.g., 100 sec—can of course, 
be performed practically with mechanical integrators 
(using a servomotor technique), but the required 
accuracy may not be economical to obtain. The first 
integration involved in going to the third-order system 
of Fig. 10 must be present anyway to obtain velocity 
from acceleration during the guidance mode. (If 
velocity meters are used this integration is inherent in 
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the instrument.) But whether to include additional 
precision mechanical integrators purely for alignment 
purposes is a problem in weight and dollar economy. 


Conclusions 


An inertial guidance system can be aligned by gyro- 
compassing, using the system’s own instruments and 
computer, and involving no external equipment (except 
for velocity information, if the base is traveling). 

Gyrocompass alignment will lead to navigation errors 
comparable with others inherent in the inertial system, 
provided sufficient alignment time and good base 
velocity information are available. 

If rapid alignment is desired in the presence of dis- 
turbing base motions, then a trade-off must be made 
between speed and accuracy. 

Fig. 10 represents a simple, practicable third-order 
gyrocompassing scheme having filtering capability 
represented by the dashed lines of Fig. 7, and producing 
speed of alignment as shown by the dashed curves in 
Fig. 9, in the absence of system nonlinearities. 

When nonlinearities are present the performance of 
the system may be deteriorated. In such cases special 
sampling techniques may produce better performance 
than the continuous system of Fig. 10. 

But, in any case, the curves of Figs. 7 and 9 represent 
an upper bound on the performance available from a 
third-order scheme. 

Improved performance may be achievable by the use 
of higher-order systems, but such systems may be quite 
expensive. 

(It should be mentioned that in specific launching 
operations it may be possible to accomplish postlaunch 
alignment, thus further reducing the required preflight 


time.) 


Appendix A 


Gyro Dynamics for Gyrocompassing Studies 


Refer to Fig. 2 which shows the stable platform, and 
to Fig. 3 which shows a unit sphere whose center is the 
center of the stable platform. Two reference planes are 
shown: a locally horizontal plane (NW) and a locally 
vertical plane (VV) which passes through the north- 
south axis of the earth. The platform is assumed at all 
times slaved exactly to its gyroscopes, so that axes 
x, y, represent both platform and gyro orientation. 


A(s) 





SKEeTcH A. (Left). SkEeTcH B. (Right) 


Y Gyro 

The spin axis of the y gyro is along x, which is located 
by Euler angles ¢, and ¢,, as shown. These angles are 
assumed to be small. 

We want to form a torque balance about the output 
axis of the y gyro, which lies in the z direction. That is, 
we want to write the relation: 


(Ay), = (My), 


in which Ay is the total angular momentum of the y 
gyro, and (My), is the resultant torque applied about its 
z (output) axis. 
The total angular velocity* of the xyz coordinate 
system—..e., of the platform—1s (from Fig. 3) 
Qp = Qe + ly ¢: + 1,4, 
in which Q, is the angular velocity of the earth: 


le = LyQz sin \ + LyQe cos A 


~*~ 


(1, is a unit vector along the y-axis, etc.) : 
For brevity, let Qg cos \X = ,, the “‘local earth rate.’ 


’ 


The components of 2p along x, y, and z are: 


lp = Q. — o,(b2 + Qe sin A) 


~ 


Qp), = 1, ” 
(Qp)y _ l, , Qp = dy — ¢, Y, 


dp = @ + Qez sin A + Py Q, 


a 
™~ 
v 
+ 
: 
II 
oe 
° 
«“ 


The total angular momentum of the y gyro rotor- 
plus-inner-gimbal is: 

Ay = 1, [Hh + I,(Qp)2) + 1, Uy(Qe),) + 12 U2(Ge)2) 
in which 7 is the spin momentum of the gyro rotor rela- 
tive to the inner gimbal, and the 1’s are the moments of 
inertia of the rotor-plus-gimbal. 

To find Hy, we apply the usual relation: 
Ay = (0H,-/dt) + Op x Ay 

(in which the partial derivative implies unchanging unit 
vectors). 

From this we obtain: 


(Ay). = T(bz + by®) + Iy(Gy — 22) [Qe — by (bz + 
Qe sind)] — [H+ 1, (Q. — dd: — oMQesin d) | (dy — o2) 
(Ay), & Iz (be + $2) — H(dy — 2) 

¢; is small, in practice, and H > JQ,, so that 
(Hy). = — H(¢, — o2,) 
was used in relation (1) in the paper. 
Z Gyro 
Similarly, for the z gyro, which has its spin axis along 
y and its output axis along x (Fig. 2): 


(Hz) _ 1, (7,(Qp).] + 1, \H so I, (Qp),) + 1, I, (Qp)z] 


so that: 


* A nontraveling base is assumed in this derivation. 


(Continued on page 912) 
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Velocity Overshoots in Magnetic Boundary 
Layers 


Paul S. Lykoudis 

Graduate School of Aeronautical Engineering, 

Cornell University, Ithaca, N.Y. 

December 5, 1960 

BOUNDARY-LAYER THEORY velocity overshoots have been 
found for the case of favorable pressure gradients and heated 

walls. Their existence has been shown mathematically by ob- 

taining asymptotic solutions at large distances from the wall 

and proving that the velocity at the edge of the boundary layer 

is approached from above.! In terms of the similarity function 


f(y), denoting the ratio of the local velocity to the velocity pre- 


vailing at the edge of the boundary layer, nondimensional ve- 
locity profiles are shown in Figs. 1 and 2. 

The question of the existence of such overshoots in the case of 
magnetic boundary layers with constant electric conductivity 
has been answered in Ref. 2. Within the frame of similarity 
solutions the conditions are exactly the same as those of the non- 
magnetic case. 

The interesting result has been found’ that for variable elec- 
trical conductivity inside a magnetic boundary layer, velocity 
overshoots are possible even in the case of wall temperatures 
cooler than the free stream. Here again the conditions for an 
overshoot were established‘ by first finding an asymptotic solu- 
tion to the differential equations of motion and energy for large 
distances from the wall. The aim of the present note is to pre- 
sent the same result through much simpler arguments and also 
offer a physical explanation. The 
It is assumed that 


same notation as in Refs. 3 
and 4 will be used. 
(1) 


= (h/he)" = gn) 


h the stagnation enthalpy, 


a/c, 
where @ is the electrical conductivity, 
n a positive number, and e denotes the edge of the boundary 


layer. The similarity equatior of motion is 
f'' +97’ = elf? — gl + QU1 — gy’)]} (2) 
where 
Qv = o-Bo?[pe(due/dx)] ~! (3) 
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Now let the subscript m denote the value calculated at the point 


where the maximum velocity occurs. The conditions there are 


i," = @ to" =e, i >a (4) 
Eq. (2) becomes 
fn’ = Blin’? + g'**O.fm' — gm(1 + Qr)] <0 (5) 


By inspection the quadratic has two roots, one positive and 
one negative, and for 8 > 0 the above inequality can only be true 
if f,’ takes a value lying between the two roots. Since f,,’ must 
also be higher than one, the positive root of the quadratic must 
be also higher than one. Hence, 


—gm'*Q,/2 + V [(gm!**Qv)/2]2 + gm(l + Qe) >1 (6) 


am 


After some manipulation the above reduces to 


&mOX gm" — 1) — (gm — 1) <0 (7) 

or 
(gu — 1)lenQOden”®™* + En" 2 +... +1) — 1) <0 @) 
For “‘cool’’ walls, gm — 1 <0 and hence the quantity in brackets 


must be positive. Since each of the ” terms inside the paren- 
theses is smaller than one, we have the two conditions 

fa® + on" * TF... ten SB (9) 
and Okan” + ga” * +... + En) > 1 (10) 


These two are compatible only if 


Qyn> 1 (11) 
For ‘“‘hot” walls with g, — 1 > 0 the reverse sign of the above is 
true: 

Quan <1 (12) 


These are the same conditions found in Refs. 3 and 4 with the 
help of an asymptotic expansion. 
From a physical point of view the above mathematical findings 


may be explained as follows. For a given value of the parameter 
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Q, (independent of the particular law of variation of o), the ve- 
locity at the edge of the boundary layer will be diminished from 
its nonmagnetic value by a factor depending only on Q,.2 On 
the other hand in the vicinity of the wall the magnetic velocity 
profile will approach the nonmagnetic one as n takes higher values 
(n > 1), since the electrical conductivity becomes smaller. For 
cool walls the enthalpy profile corresponding to the case of vari- 
able electric conductivity will lie above the one calculated for a 
constant value (Fig. 1). In this case the mass density is lowered, 
and hence flow accelerations resulting in velocities above the ones 
prevailing at the edge of the magnetic boundary layer may occur 
The higher the power m in Eq. (1), the larger is the distance over 
which the velocity profile will be closer to the nonrnagnetic one, 
and overshoots may be observed above smaller values of the 
parameter Q,, consistent with the condition of inequality [Eq 
(11)]. 

In the case of hot walls the situation is different. 
it can be seen that the effect of an increased magnetic field is to 
increase the enthalpy profiles (same effect as in the case of ad- 
Assuming that we start from a situa- 


In Fig. (2) 


verse pressure gradients). 
tion in which a velocity overshoot is possible without the presence 
of a magnetic field, such an overshoot will be conserved in the 
magnetic case if the electrical conductivity remains constant 
For variable conductivity, however, the enthalpy profile will fall 
below the one calculated with constant conductivity, the mass 
density will be higher, and thus local decelerations will occur. 
It is also evident that the original velocity overshoot will disap- 
pear above a critical value of Q, which will be smaller for larger 
values of n, a fact consistent with inequality [Eq. (12)]. 
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Simple Approximate Solutions for Tangential, 
Low-Thrust Trajectories 


Harold M. Stark and Paul D. Arthur 
Systems Corporation of America, Los Angeles, Calif 
January 16, 1961 


SYMBOLS 


R = distance from origin 
V = velocity 
6 = central angle (true anomaly) 
y = inclination angle with respect to local horizontal 
g = acceleration due to gravity 
t = time measured from the initial conditions 
a = constant acceleration, made nondimensional by go 
7 = nondimensional time, 
3 Vol sin yr 3Voat 
= 1 + — = ] a = [Eq. (7)] 
2 Ro Ro 


Subscripts 
r = reference orbit 
0 = initial conditions 


I OW-THRUST TRAJECTORIES are of interest for interplanetary 
travel; however, exact solutions for these trajectories do not 


in general exist. Thus there is a need for simple approximate 
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solutions; they have been found for constant combinations of 
radial and circumferential thrust and have been verified for the 
case of constant radial thrust.!. For the case of constant tangen- 
tial thrust, trajectory characterizations are presented here which 
are shown to be valid up to the last few revolutions before escape 


DERIVATION OF EQUATIONS 

The usual inverse-square central gravitational force field is 
assumed. In a previous work! deviations from a reference 
Keplerian orbit due to application of small constant radial and 
circumferential thrust were calculated. Here, however, due to 
the spiral nature of a low-thrust, tangential trajectory, deviations 
are taken from a logarithmic spiral 6, ~ In (R,/R), using a non- 
constant tangential acceleration of ag. For easy reference, sev- 
eral pertinent equations for logarithmic spirals are presented: 








V; = —1/2g, sin yr (1) 
Vevr — (Vr?2/Rr) COS yr = —£r COS Yr (2) 
<i 0 (3) 
R,/Ro = &e0 vw = V3 (4 
V; grR, = Vor—"3 (5) 
R,76, = VoRor? cos yr (6) 
sol 
—— EXACT 
——— £Q (14) 
r / 
ae / 
Ro | 
RE + ° 
20 es, Fic. 1. Approximate and 
exact trajectory characteris- 
- . ef 
tics for a = 0.001. 
| 4 
L 
r 5 REV 
1o 4 eS 1 
0 100 200 300 400 
Vot 
Ro 


sin yr = 2a (7) 


where dots indicate differentiation with respect to time. Eggs. 
(1) and (2) are the equations of motion of a trajectory in tangen- 
tial and normal coordinates with a tangential thrust equal to 


a 
58r Sin Ye and zero normal thrust. Egs. (3) through (6) give 


the characteristics of the set of logarithmic spirals satisfying Eqs. 
(1) and (2). In Eq. (7) that particular value of y, is chosen 
which gives the reference orbit an initial acceleration of ago. 

The equations of motion in radial and circumferential com- 
ponents of the actual and reference trajectories are 


R — R(6)? = —g + ag sin y (8a) 
Re — R,-(6,)? = —gr + ag sin yy (8b) 
(1/R)[d(R%)/dt] = ago cos y (9a) 
(1/R,)[d(R,%,)/dt] = ag, cos vr (9b) 


For low-thrust trajectories, a < 1. This implies y, < 1. In 
addition assume that y and the quantity x, defined by R = R,(1 
+ x), are small. Henceforth, products of small quantities will 
be neglected. The only exception is aV¢t/Ry which is less than 
one in the range of ¢ considered here but is not as small as x or a. 
Thus (aVot/Ro)? will be kept but higher-order terms will be 
neglected. Utilizing these assumptions, Eq. (8b) is subtracted 
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from Eq. (8a) to give 
d?(R,x) (R%)2(1 — 3x) — (R,6,)?  2goRo? 
= - ea =i (10) 
dt? R,3 R,? 


For constant a, Eqs. (9a) and (9b) can be integrated and ex- 
panded: 
R%™ — R,6, = 2VoRo(aVot/Ro)? (11) 
Using Eggs. (11) and (6) in Eq. (10) results in 
d?(R,x)/dt? = (goRo?/R,?)[—x + 4(aVot/R 24/R R,] (12) 
The solution to Eq. (12) is 
x = (4/9)r-"38(7 — 1)? — 2ar- sin 6, (13) 
R/Ro = 7r¥%* + (4/9)r"3(7 — 1)? — Zar’? sin 6, (14) 
The velocity V can be found from an energy balance. The 
final energy of the orbit equals the initial energy of the orbit 
plus the work done by the thrust: 


V2 = goRo? goRo : “ 
2 ‘ R aplee 2 i ago Vat’ (15a) 
7 V 
V,? goR goRo ; = 
— sms = —- r ag, V-dt' (15b) 
2 R; R, 6 
Now if we define y by V = V,(1 + y) and assume y < 1, then 


subtraction of Eq. (15b) from (15a) yields 


Vi2y = goR (® =) + : = 6 
fy = goRo | > R, ago : V.11- E) dt' (16) 


Integrating via 7 and using 1/(1 + x 1 — x gives 
l i. : = 
yso—xtorls — 723 (17) 


The assumption that y and y, are small results in 
6= V/R, 6, = V,/R, (18) 
The difference of 6 and 6, can be expressed as 


6 — 6, = (Vo/R ry — x (19) 


Integrating Eq. (19) and adding 8, gives 
@ = 1/a [(1/2) In r — (49/40) + (8/9)r-"3 +(7/18) r 78 + 
(1/8)r#’3 — (8/45)r*/3] + 4a (1 — 77“ cos 6,) (20) 


In all cases where sin 6, and cos 6, appear, the terms involved 
are unimportant for large ¢ but are often the dominating quanti- 
ties for small ¢. For small ¢, the quantities V and @ reduce to the 
expected: 

V = Vo + agot; 6 = (Vot/Ro) + (aget?/2Ro) (21) 


These equations are in complete agreement with the circumfer- 
ential equations of Ref. 1. This is of course expected since ini- 
tially tangential and circumferential thrust trajectories are 


identical. 


COMPARISON OF APPROXIMATE AND EXACT TRAJECTORY 
CHARACTERISTICS 

Analytic solutions for a constant-tangential-thrust trajectory 
cannot be found. Hence computer solutions? must be used for 
comparison purposes with the approximate characteristics de- 
veloped here. Exact and approximate solutions of R as a func- 
tion of ¢ for a = 0.001 are presented in Fig. 1. It can be seen 
that the approximations are accurate for about thirty revolutions. 
For a = 0.001, escape occurs at 39.95 revolutions. Approxima- 
tions for V and @ have produced similarly accurate results. 
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A Simplified Method for Computing 
Real-Gas Inlet Efficiencies 


C. W. Hartsell* 
The Aerospace Corporation, Los Angeles, Calif. 


January 16, 1961 


SYMBOLS 


1 = integration constant 

Cp = specific heat at constant pressure 
h = enthalpy 

KE = kinetic energy 

VW = Mach number 

P = pressure 

R = gas constant 

5 = entropy 

7 = temperature 

if = velocity 


n = kinetic-energy efficiency 
= ratio of specific heats 


Subscripts 
= free-stream ambient conditions 
= conditions at the end of the compression process 


2 = 
s = condition at Po upon isentropic expansion from station (2) 
T = stagnation conditions 

e = isentropic path 


Bb PERFORMANCE of inlets at hypersonic flight speeds is com 
plicated by the required use of numerical techniques. A 
basic computation is the evaluation of the inlet kinetic efficiency 
(yn). The evaluation of y for real-gas flows is most easily described 
with the aid of the Mollier-diagram schematic in Fig. 1. The 
inlet is operating with free-stream ambient conditions signified 
For an adiabatic system the inlet compression process 
) —> (2) path where Point (2) is at the free- 
To evaluate the inlet, theoretically 
Point (3) 


as (.). 
occurs along the 
stream total enthalpy (/7 

re-expand the gases isentropically from (2) — (3) 


is determined as being at the free-stream ambient pressure. The 
kinetic efficiency is now defined as follows: 
0 (hr — hs)/(hro — ho), = KE2:/KE (1) 


Thus, the ratio of the kinetic energy of the captured gases when 
expanded to the ambient pressure divided by the free-stream 
kinetic energy defines the system efficiency. 

A number of practical problems arise from evaluation of 7 as 
For high-efficiency inlets 7 is close to 1 so that 


described above. 
The evaluation of 


high computational accuracy is required. 
Point (3) is usually difficult as it is usually found in a low-accuracy 
portion of the Mollier diagram. A little experience with this 
technique for a large number of practical cases discloses that 
Point (3) is often located in a region of no dissociation and where 
yis close to7/5. With these restrictions it is possible to combine 
the first and second laws of thermodynamics to yield an expres- 
sion for the entropy rise along the constant pressure line (P.). 
Thus we have 

S = Cpln T + a (perfect gas) (2) 
then the entropy increase between ( ~ ) and (3) results in 

AS = Cp In (7;3/T.) ( 

Taking the definition of 7 we may insert Eq. (3) as follows: 


hre — he + he — h 14 CpeT.(1 — T3/Ta) (4) 


——© 


_ KE KE., 
then 
sai Cri al St sated (5) 
- AS 7-1 
e= 1 = Bly ~~ 2DeM CE Ye =» (6) 
if + = 7/5 


* This work was initiated while the author was connected with 





The 


Marquardt Corporation. 
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where AS Ss —- § 


A direct check of Eq. (7) for the normal-shock case against the 
exact solution is presented in Fig. 2. Ref. 1 presents the results 
of normal-shock real-gas computations which define the Point (2) 
conditions. The Mollier data for air are found in Refs. 2 and 3 
It will be noted that this approximate solution yields correct re 
sults to a Mach number of 10 for both altitudes studied. It 
would be expected that for a typical high efficiency inlet (n > 
0.90) the region of usefulness of this approximation would be ex- 
tended considerably above Mach 10 

The usefulness of the above technique is expected to be in the 
region of moderate real-gas effects. For the case investigated 
herein this region is presumed to be in the range of 5 << M, <10 
where the inlet efficiency is within 10 percent of the perfect-gas 
It is apparent that most present work lies in this flight- 


value 
In this region the present technique provides rapid 


speed range 
and consistent answers with very little if any reduction in ac- 
curacy. In contrast, the ‘‘exact’’ technique has been found to 
yield reasonably consistent answers only after very careful and 


tedious cross-plotting. 
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Fic. 2. Evaluation of approximate solution (normal shock 
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Thermal Stresses in an Idealized Wing 
Structure 


Mark Levinson 

Assistant Professor of Mechanical Engineering, 

Oregon State College, Corvallis, Ore. 
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vers THERMAL STRESSES due to aerodynamic heating 
of the idealized wing structure shown in Fig. 1 have 

been computed, by means of an approximate two-dimensional 
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Fic. 1. Idealized wing structure. 


heat conduction analysis, for the particular case called Model 5 
be Tramposch and Gerard.! The heat-conduction analysis is 
an extension of Biot’s? work on the one-dimensional problem. 
Since the problem is of rather academic interest, the main pur- 
pose of this note is to indicate the usefulness of Biot’s variational 
principle of heat flow. No use is made of classical results. 

Surface temperature of the isolated flange—i.e., in the ab- 
sence of the web—is assumed to be 

T; = To[1 — exp (—8t)] (1) 
where the constant 8 is found from the law of conservation of 
energy and 7) is the adiabatic wall temperature. Thus B = 
K/ac where K is the boundary-layer heat-transfer coefficient and 
c is the heat capacity per unit volume of the structure. It is 
possible now to compute temperature distribution in the flange 
for all time. Only times greater than the transit time in the 
flange, /;, are of interest 

Tis(x, t) = (Ty — To) [1 — (x/a)]}? + Tr (2) 
where 
Ty = To{1 — exp[n(t; — t)] + (4/17)8 + n/B — n exp X 
(—Bt) [exp {B(4 — t)} — exp {n(ti — t)} ]} (3) 
and n = 42/17 (x/a?) where x is thermal diffusivity. 

For ¢ > ¢, there is an interaction between the web and the 
flange. The effect of the web is to act as a heat sink and to cool 
the portion of the flange near the junction of the two elements. 
Temperature distribution in the affected portion of the flange is 
given by 

Tas = Tis{1 + (Mm — 1)(¥?/Q?)] (4) 
where q is a generalized thermal coordinate in Biot’s sense, and 
q is the penetration depth given by Biot as 

q = 3.36[(t — t)]!/? (5) 
Tay is shown graphically in Fig. 2. 

Temperature distribution in the web is given, until tempera- 
ture starts to rise at the middle of the web, by 7, = 9: 7%(%1?/q?). 

By means of Biot’s principle and the law of conservation of 





Fic. 2. Temperature distribution in the portion of the flange 
affected by the web. 
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energy, a first-order, linear, ordinary differential equation with 
quite complicated time-dependent coefficients is obtained for 
gq. The initial condition is that gi(4;) = 1. For the particular 
case presented in this note, the equation was solved numerically 
with the ALWAC 3-E computer owned by the Department of 
Mathematics at Oregon State College. 

After a time ¢ = t; + te, where tf, = 0.0885 1,?/x, the tempera- 
ture at the middle of the web starts to rise. The simplest way 
to determine the temperature at the middle of the web after it 
starts to rise is to apply Duhamel’s integral to Biot’s solution of 
the problem of a slab with constant temperature on one surface 
and insulated on the other surface. 


“4 —0.214 1(qi:Ty) 
IT. = f J; — exp (t'’ — r) ( = dr (6) 
o ) to f dr 


for t = bh, + ty + fe. 
Temperature distribution throughout the web for ¢ > 4; + tis 
given by 
Tw = (M75 — Tm)(%1?/h?) + Tm (7) 
It is now possible to calculate the thermal stress at the middle 
of the web as a function of time by means of the equation 
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Fic. 3. Thermal stress at the middle of the web (model 5). 


o = —aETy + aay f, aETdA (8) 


where the integration is performed over the cross-sectional area 
and a@ is the coefficient of linear thermal expansion. Explicit 
mention should be made of the fact that a temperature-dependent 
modulus of elasticity is used in the numerical example considered 
in this note. 

The physical and geometrical parameters for the case presented 
in Fig. 3 are those of Model 5 in Ref. 1. Similar results were 
obtained for Model 6. Comparisons are made with the experi- 
mental and theoretical results given by Tramposch and Gerard 
as well as the earlier one-dimensional theories of Torda and Hoff’ 
and Pohle and Oliver. 

Since the boundary-layer heat-transfer coefficient measured 
by Tramposch and Gerard varied considerably during the 
course of the experiments, and all the theories assume the co- 
efficient to be constant, it is necessary to comment upon the effect 
of this variation. The maximum value, at the start of the ex- 
periments, was about 19 and the final value was about 11. It 
seems reasonable to believe that this would result in the thermal 
stresses rising more rapidly at first and then decreasing more 
rapidly than if a constant heat-transfer coefficient had been 
maintained. If this is so, it constitutes an explanation of the 
somewhat different shape shown by the _ stress-versus-time 
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curve predicted by the present theory and the one found experi- 


mentally. 
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Hot-Wire Measurements in the Hypersonic 
Wake of a Cylinder} 
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¢ rs WAKE behind a body moving at hypersonic Mach num- 
bers is an interesting and largely unexplored type of free- 
shear flow. In order to acquire an understanding of the general 
features of such a flow, the wake of a circular cylinder positioned 
normal to a flow of air at a Mach number of 5.8 is under investiga- 
tion in the GALCIT hypersonic wind tunnel. This note con- 
tains preliminary information obtained with a hot-wire ane- 
mometer used in this wake 

Perhaps the two most immediate questions arising in this 
problem are: (1) the conditions under which the wake is laminar 
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Fic. 1. Schematic diagram and nomenclature of the two 
dimensional wake behind a circular cylinder. 


~ 


region of shock-heated gas trailing behind the cylinder. Both 
of these questions can be answered with reasonable precision 
by using the qualitative diagnostic properties of the hot-wire 
anemometer, without resolving the data into temperature or 
mean-flow fluctuations and their spectral distributions. For 
example, by moving a constant-current hot wire from a region 


or turbulent, and (2) the rate of growth of the wake into the 


Tt The work discussed in this paper was carried out under the sponsorship 
and with the financial support of the Office of the Chief of Ordnance, Army 
Ballistic Missile Agency, and the Office of Ordnance Research, U.S. Army, 
Contract No. DA-04-495-Ord-1960. 

* Senior Research Fellow. 
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Fic.2. Effect of the free-stream Reynolds number per inch Re’ 
on the position of the transition point on the center-plane of the 
wake at M = 5.8. The ordinate of each plot is proportional to 
the mean-square output of the wire; the maximum on each plot 
marks the transition point. 


of laminar to one of turbulent flow one observes a characteristi- 
cally large increase in the mean-square wire output. Even though 
exact statements about the turbulent intensity cannot be made 
in such cases unless corrections for the wire sensitivity changes 
and thermal lag are made, one can at least by comparison dis- 
tinguish turbulent boundaries (such as boundary-layer or 
wake edges) and transition regions 

Measurements based on this principle were made at fixed 
wind-tunnel supply pressure by moving a constant-current hot 
wire in a streamwise direction along the wake centerline [the 
mid-span line of the plane defined by the cylinder axis and the 
flow direction (Fig. 1)]. Typically, the mean-square output of 
the hot wire increased for some distance beyond the cylinder, 
attained a maximum, and then decreased monotonically over 
large distances downstream, as shown in Fig. 2. Such behavior 
is explainable as a true qualitative variation of the turbulent 
intensity. Consequently, the maxima of these streamwise 
traverses were defined as the ‘‘transition’’ points (represented 
by X,.). The pretransition signal increase is interpreted as 
the gradual increase of flow intermittency which precedes 
transition, and the posttransition decrease in the signal is inter- 
preted as turbulence decay. 

The transition data obtained with cylinders 0.300 in., 0.200 
in., and 0.100 in. in diameter (corresponding to cylinder Reynolds 
numbers ranging from 15,000 to 56,000) are plotted in Fig. 3, 
where Re’ represents the Reynolds number per inch in the free 
stream preceding the bow shock The most striking result is 
the absence of any effect of the cylinder diameter. In other 
words transition here exhibits a ‘‘flat-plate’’ effect, where the 
dimensional transition length depends on the ambient conditions 
only, and where the history of the wake up to the formation of 
the ‘‘neck”’ has little effect on the transition point, at least within 
the range of cylinder Reynolds numbers investigated. The 
transition Reynolds number based on the free-stream conditions 
(and the transition length) was found to be of the order of 
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Fic. 4. Growth of the nondimensional wake thickness b/D with 
the nondimensionalized streamwise coordinate x/D (Fig. 1) for 
various cylinder Reynolds numbers Rp. 


300,000, while the Reynolds number based on quantities at the 
boundary of the wake at the transition point was of the order 
of 60,000. 

Another measurement was made aimed at studying the 
variation of the wake thickness 6 with x. For this purpose the 
wire was traversed along y (Fig. 1) for a number of x’s and tunnel 
supply pressures. From the observed variations of the mean- 
square output with y it was possible to discern the wake bound- 
aries and thus define an approximate wake thickness. This 
procedure was possible even in the case of laminar wakes, where 
a small amount of intermittency in the wake resulted in measure- 
able rms signals, but not so for the smallest cylinder Reynolds 
numbers (below about 6,000). It is not clear, at this time, 
whether this latter effect is connected with the stability of the 
wake or is a limitation of the experimental equipment. 

Fig. 4 shows the results obtained when the nondimensional 
wake thickness })/D is plotted as a function of x/D. The data 
cover a range of cylinder Reynolds numbers from 6,500 to 71,000, 
but show practically no dependence on that parameter. In 
fact, the points seem to fall generally along two lines, correspond- 
ing to laminar and turbulent wakes (cf. Fig. 3). In both cases 
the thickness }/D appears to grow as some power of x/D ap- 
proximately equal to '/s, at least up to x/D of about 100. These 
effects are now being studied in more detail. 


Additional Information on the Calculation of 
Natural Modes of Free-Free Structures 


James A. Stricklin 
Assistant Professor of Aeronautical Engineering, 
Georgia Institute of Technology, Atlanta 13, Ga. 


February 10, 1961 


SYMBOLS 
C(x,9; &,n) = influence function 
k(x,y; &n) = stiffness function computed with three points attached 
&1(t) = generalized coordinate 
p(x,y) = mass density per unit area 
w = natural frequency 


(other symbols are the same as in Ref. 1) 


_— has arrived at the interesting result! that in calcu- 
lating the free-free modes of a structure one may use in- 
fluence coefficients for the body simply supported at three 
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points—or fixed at a single point—and rigid-body modes with 
respect to the center of gravity. It is the purpose of this note 
to give further insight into the method through the use of assumed 
modes and Lagrange’s equations. 

The discussion will be restricted to a plate-like structure 
simply supported at three points whose only deflections of 
interest are those in the z direction. 

Assume a set of functions f;(x,y), i = 1, 2,..., 
the property that any deflection consistent with zero deflection 
at the support points an be represented by a linear sum of these 
functions. Each of these functions must correspond to a possible 
displacement of the structure simply supported at the three 
points, but need not satisfy the nongeometric boundary condi- 


©, which have 


tions.* 

Then any deflection of the structure can be represented as a 
linear superposition of the f,;(x,y)’s and three additional linearly 
independent functions which free the attach points. The three 
additional functions may be entirely general’ in that they can 
include elastic deflections of the structure; but this would only 
create additional work for the analyst. Therefore, let the three 
functions be a rigid-body translation and two rigid-body rota- 
tions about principal axes through the center of gravity. 

Thus, any deflection can be represented as 


w(x, y, t) = wo(t) + FOr(t) + £0,(t) + > filx, y)é&(t) (1) 
$=] 


The internal energy and kinetic energy are given by 


1 2 e x 
v= ff, fe ye {fe 98" x 
‘= i=] 


Filé, n)&i(t)dtdndxdy (2) 
l oe 
i ff, [(x, y, t)]* p(x, y) dxdy (3) 


Applying Lagrange’s equation (which is true for any virtual 
displacement**) to the generalized coordinate £;(t) and manipu- 
lating OU /O£; by interchanging order of integration and dummy 
variables of integration, remembering that k(x,y; §,n) = k(é,n); 
(x,y), one obtains 


ff, ‘ff, k(x,9; §,n) a SilEn )E;(t) didn + 


w(x, y, t)p(x, yt Ix, ydxdy = 0 (7 = 1,2,..., @) (4) 


Then the only way this can be satisfied (except at an attach 


point) for all 7 is for 


ff, k(x, y; &, n) ys F.C, néi(t)didn + w(x, vy, t) p(x,y) = 0 
t=1 
(5) 
or 
SSe k(x, y; &n) > SilEm)E(t)dtdn = — w(x, y. t)p(x, y) (6) 
con 
where S’ excludes the attach points. Eq. (6) follows from Eq. 
(5), since it is possible to exclude the attach points from the 


original expressions for kinetic and potential energy. Inverting 
Eq. (6) and assuming w(x, y, ¢) to be simply harmonic one gets 


w(x,y) — Wo — XOy — FOr = w? f fice: y; &n) p(é,n)w(é,n )dédn 
(7) 


Applying Lagrange’s equations to the rigid-body degrees of 
freedom, one obtains the equations for conservation of linear and 


angular momentum: 


ff p(x, y)w(x, y, t)dxdy = constant = 0 (8) 
Sf, p(x, y) w(x, y, t)dxdy = constant = 0 (9) 


ff, p(x, y)pw(x, y, t)dxdy = constant = 0 (10) 
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{(P) = F'(P) T cos y cos @ (2) 
ual I ET THE WING be constructed such that straight generators The total strain energy V of the system takes the form 
yu- form its surface. The surface is a cone with its base repre- . 
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n = [(T cos y cos ¢ sin? 0)/Ipx](QX?/x) (4) 
where Jpx is a reduced second moment of area of the root sec- 


tion 
L , ms 
Igy = ‘f, T? cos? y cos? ¢ sin? 6 dP (5) 


Eq. (4) can be written in the form 
n = [(x/X)T cos y cos ¢ sin? 6)M/Irz (6) 
where Jp; is a reduced second moment of area of the section at 
distance x 
Tre = (x/X)Trz (7) 
As an example, calculations are performed on the delta wing 
shown in Fig. 2; the values of ”/nmax are also plotted. 


CONCLUDING REMARKS 
A question now arises whether Eq. (6), together with Eq. (7), 
could represent, in the case of sweptback and delta wings, the 
version of the bending equation of simple beams. The latter equa- 
tion can be deduced from the former in case y = ¢ = O and 6 
= 90°. The author suggests an affirmative answer. 
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Method for Calculating Surface Streamlines 
and Laminar Heat Transfer to Blunted Cones 
at Angle of Attack} 


E. A. Sanlorenzo* 
General Applied Science Laboratories, Inc., 
Merrick and Stewart Avenues, Westbury, N.Y. 


February 20, 1961 


NOMENCLATURE 


p = pressure 

AS/R = nondimensional entropy 

Y = ratio of specific heats 

Vr = radial velocity component 

Ww = cross-flow vi locity component 

0 = polar angle about the longitudinal axis of cone measured from 
the most windward meridian 

V = velocity, nondimensional with respect to the limiting velocity 

a = angle of attack 

s = distance along surface measured from a = 0° stagnation point 

Ro = radius of spherical nose 

Subscripts 
se = stagnation at edge of boundary layer 
© = freestream 


AGLIO-LAURIN! provides a method for calculating the 
laminar heat transfer to general three-dimensional blunt- 
nosed bodies. Following this approach, after a suggestion of 
Dr. A. Ferri, a relatively simple method has been obtained which 
allows the determination of the heat transfer to a blunted cone 
at angle of attack if the pressure distributions along three me- 
ridionalrays are known. The method, including a more extensive 
comparison with experiment than is possible here, is given in 
Ref. 2, the principal results of which are summarized here. 
In Refs. 3 and 4 it was shown that the peripheral pressure dis- 
tribution for the cone at angle of attack, a, can be closely ap- 
proximated by a cosine distribution of the type 


b/pse = (b/Pse)aeo + Aa cos 0 + Ba? + Ca? cos 20 (1) 


where the empirical coefficients A(s), B(s), C(s) are determined 
from the pressure at three peripheral locations at each station 
along the longitudinal axis of the cone. 

To determine the heat transfer at any point on the cone sur- 
face, the inviscid streamlines over the body must first be deter- 


{ This work was sponsored by WADD, Flight Dynamic Laboratory, 
Contract No. AF 33(616)-6692, and is part of a larger theoretical and experi- 
mental program. The experimental data were obtained at AEDC, Tulla- 
homa, Tenn. 

* Senior Scientist. 
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mined. The heat transfer is then calculated by an integration 

along the streamline taking into account the local fluid properties 

and the spreading of the flow about the particular streamline. 
The location of the streamlines may be approximated by as 


suming a polar coordinate system with the velocity components 


on the surface given by 
V, = Vo + Vi cos 0 


(2 
W = W, sine f 


where Vo, Vi, W; are coefficients which must be determined at 
each peripheral station of the cone. 

At any point on the cone surface the pressure and velocity are 
related by 

b/po = [(1 — V2)/(1 — Va2)]¥/7~! e—45/R (3) 
For a blunted cone having a moderate ratio of base radius to 
blunting radius (<10), the entropy at the edge of the boundary 
layer may be taken to be the entropy rise through a normal 
shock. In this case, if the entropy rise is assumed constant, 
and the pressure along three meridional rays of the cone is 
known, Eqs. (1), (2), and (3) may be combined and values of 
Vo, Vi, W; obtained at each peripheral station of the cone. 

At any point of the body (corresponding to a particular sta- 
tion and a particular value of @), it is possible to compute V, and 
W using the set of coefficients Vo, Vi, W; which are appropriate 
to that station. Since W is normal to V;, the direction of the 
streamline at a point is given by 

e = tan! W/V, (4 

By proceeding in a stepwise manner, a previously computed 
value of e« permits the determination of the streamline slope and 
the location of an adjacent point at the next chosen station, 
where the computation is repeated using the value of @ and the 
set of coefficients Vo, Vi, W; appropriate to the second station. 

Using the notation of Ref. 5, the heat-transfer parameter 
Nu/V Re is given by 

Nu fart’ *..__. l a ee 
—=> = 0.353 (Pr)!/3 {| —— (piitthz) ppitthe?dl 
V Re Pre 0 


(5) 


Eq. (5) must be integrated along the streamline with the quan- 
tity pia determined from the pressure on the body and the stag- 
nation pressure at the edge of the boundary layer. From the 
streamline pattern and Eq. (1), the pressure along any stream- 
line can be calculated. The quantity /, which characterizes 
the spreading of the streamline at a particular point is also deter- 
mined from the streamline pattern by summing the length 
elements between the point on a streamline and two adjacent 
streamlines. The length elements are measured from any point 
on the streamline in a direction perpendicular to the adjacent 
streamline. The adjacent streamlines are chosen so that the 
mass flows contained between the particular streamline and the 
adjacent streamlines are equal. Proceeding from point to point 
along the streamline, the variation of the spreading can be deter- 
mined. To insure compatability with the heat-transfer results 
over the sphere, the variation of 42 along the streamline must be 
scaled using the spreading at the sphere-cone junction corre- 
sponding to the flow over the sphere. At the sphere-cone junc- 
tion the streamline spreading is proportional to the sine of the 
angle through which the flow turns in passing over the sphere. 

The pressure and heat transfer distribution over the sphere 
are well known—e.g., the Fay-Riddell stagnation-point theory 
and Lees’ distribution of pressure and ratio of local to stagnation- 
point heat transfer—consequently at the sphere-cone junction 
the heat transfer and other fluid properties are known and there- 
fore the integral quantity of Eq. (5) is determinable. Using 
this initial integral value the heat transfer along any streamline 
is obtained by numerically integrating Eq. (5). 

The method of analysis has been compared with experimental 
results obtained at M = 8. Details of the body geometry and 
test conditions are indicated in Figs. 2 and 3. The results of 
two cases are presented here, a = 15° and a = 30°; in both cases 
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normal-shock stagnation conditions were assumed to prevail 


outside the boundary layer. Fig. 1 shows the streamlines on 


a development of the conical portion of the 15° half-angle cone 
The heat-transfer results along 0 0), 


2 and 3 for the 


for the case of a 30° 
ind 2/2 are compared with experiments in Figs 
The experimental data were obtained over 


both 


two angles of attack 
the Reynolds number range 0.38—2.3 million per ft. In 
cases, good agreement with the experimental results was ob 


tained 
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Pressure Distribution on Slender Bodies and 
Wings Oscillating at High Frequency 


Meir Hanin 
Technion—lsrael Institute of Technology, Haifa, Israel 
March 6, 1961 


= A RECENT NOTE! an asymptotic representation was given for 
the unsteady flow produced by a slender configuration (wing, 
body, or wing-body combination) oscillating harmonically at 
high frequency. The results included the second high-frequency 
approximation to the oscillatory pressure on the surface of the 
configuration. In the following, the third high-frequency ap 
proximation to the pressure on the surface will be obtained 

As in Ref 
pendicular to the direction of flight, and that the oscillatory 
We de 


note by (, 7, ¢) a body-fixed Cartesian coordinate system in which 


1, we assume that the oscillating motion is per- 
velocity is small in comparison with the flight speed 


the ¢ axis is taken in the direction of the oscillation, the £ axis 
in the mean downstream direction, and all the distances are re- 
ferred to a standard length /] (body diameter or wing span). 
Since the configuration is aerodynamically slender, the oscillatory 
velocity potential ye'®! will obey the two-dimensional Helmholtz 
equation in the (y, ¢) coordinates [Eq. (1) of Ref. 1], and the 


iscillatory pressure pe'®’ will be related to the potential through 
= — py Iolr ly (1) 


Pp = —ipowy 


where ap and p) denote the speed of sound and the density of 
the undisturbed air, w denotes the angular frequency of the os 


cillation, and A the reduced frequency defined by 
A wl /ay 


We assume that the reduced frequency A is large. From Eqs 
(3), (4 1 it follows that the asymptotic expan 


sion of the potential for large \ has the form 


, and (10) of Ref 


y¥ _ |r 1Fa 4- X 2Fe 1. > F 4 je thy (9 


where the coefficients F\") depend on (£, n, ¢) only, and »v is the 
distance (taken along the normal in the & 
n, ¢) to the surface of the configuration 


const. plane) from 


the point (&, Having 


found already the function F“ in the flow field and the values of 


F® on the surface [cf. Eqs. (13) and (14) of Ref. 1], we can 
evaluate F® on the surface without carrying out additional 
integrations. In fact, Eqs. (8) and (10) of Ref. 1 give 

( F®), —i(O0F/dyv), ) (3 


while from Eqs. (6b), (11), and (12) of Ref. 1 it follows that 


(= . ) 1 (= : oO? FU ) l ‘ 
_ t — ( Fe), (4 
Ov /.=0 2 On? or? J, 2R 


R denoting the local (dimensionless) radius of curvature of the 


surface in the — = const. plane 
In order to evaluate the Laplacian of FC 


the 


we have to pass trom 


the (ny, ¢) coordinates to “characteristic’’ coordinates in 


terms of which F“) was determined. The expression for the 
Laplacian turns out to be somewhat lengthy, but on the surface 


it reduces to 


0? FU 0? Ft ] , d? : 
4 - (Fa : +- (FO) (5 
On? ot? Jy»  4R? 1s? 
where s denotes length (referred to /) along the ¢é const 
contour of the surface. 
The third large-\ approximation to the pressure on the surface 


can now be obtained at once by combining Eqs. (1) to (5) and 


substituting the values of PF and F® at v 0 from the results 
of Ref. 1. It is found that the unsteady pressure on the surface 


of a slender configuration oscillating at high frequency is given by 


p — polo [w cos 8B + (1/2)! Rw cos B — A7* X 
\(3/8)R~*w cos B + (1/2) [d2(w cos B)]/ds?} + O(A~§)] (6) 


Here 8 denotes the angle between the local normal to the surface 


t . 
we" is the 


and the upward direction of the oscillation, and wv 
local oscillatory velocity on the surface, w being counted positive 
when the motion at ¢ Note that d8/ds 

R-, and that —w cos8 e'®’ is the outward-normal component 


of the oscillatory velocity of the surface 


0 is downward 


In the particular case of a body of circular cross section, 
oscillating so that w does not vary along the circumference, Eq 
(6) gives 
p = —podow cos B[1 + (7/2)(AR) 4 


(1/8)(AR)~? + O(A78)] (7) 


For this case the two-dimensional Helmholtz equation can be 
solved exactly, aad we find that Eq. (7) agrees with the asymp- 


totic expansion of the exact solution 
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For a flat wing or wing panel oscillating normally to its plane The general airfoil configuration to be used in this analysis is 
we get from Eq. (6) the symmetrical, double-wedge airfoil pictured in Fig. i. Em- 
Pp = F potyl[w — (1/2)\~*(d2w/ds?) + O(A-4) | (8) ploying a two-dimensional, two-degree-of-freedom analysis, 


the equations of motion become 
where the signs * refer to the upper and the lower surface 


respectively, and s is taken along the span. mh + moy*h L 0 
REFERENCE Fate F Vegiter*er =e - 

' Hanin, M., High-Frequency Oscillation of Slender Bodies and Wings where LZ and JV are the oscillating lift and moment terms, #7 
Journal of the Aero/Space Sciences, Vol. 27, No. 7, pp. 544-545, July 1960" airfoil mass per unit span, /, mass moment of inertia per unit 
span about elastic axis, w, = natural frequency in bending, and 
+ Wa = natural frequency in torsion. Structural damping has 
been neglected and the center of gravity assumed coincident 
with the elastic axis at the center of the airfoil. Use is made of 
A Note on the Effect of a Time-Varying Forward the second-order, linearized, piston theory, as formulated by 
Flight Velocity on the Bending-Torsion Ashley and Zartarian,! in developing the expressions for the os 
Stability of a Supersonic Wing cillating lift and moment terms. Omitting the details of the 


derivations, 1 and M/ can be written as 


Robert L. Swaim, 1/Lt. USAF 





- L -~pch|4b + (y + 1)(U/c)*%r?| 4 
Flight Control Laboratory, Wright Air Development Division, ee pee jp etal oe 
Wright-Patterson AFB, Ohio peal(y + 1){U/c)b?7] — pcUal4h + (y + 1)(U/c)%br?] (8 
‘ 
February 15, 1961 M = pch\(y + 1)(U/c)b*z] — 
pea} ( 4/3)b3 + [( y + 1)/3](U/c)2b3r?} 4 
[ ‘ig TO THE PRESENT time flutter analyses have generally been pcUal(y + 1)(U/c)b?2r] } 
carried out under the basic assumption of constant flight ; ' ; : : 
i ka aie where p free-stream density, ¢ speed of sound, / free 
velocity through the critical range. This has been shown to , stabi : ; 
stream velocity, 7 maximum thickness-to-chord ratio, and 


be a reasonably valid assumption when applied to the flutter his ate ws 
: : é fH : 7 ratio of specific heats. 

analysis of vehicles whose maximum rate of acceleration through i ; : ; ’ : 

of ; ; Employing Eqs. (3) and (4), Eqs. (1) and (2) can now be written 

the critical speed range is less than approximately lg. But : 

as 


when higher accelerations are involved, the flutter speeds ob 
tained through analysis under this assumption become increas- A+ te + alk + oe + ck tte + ala 0 
ingly in error from the actual values. It was of interest to the 


author to investigate the degree of this error for the case of a a + (7 + esU*)e& + cQUh + (CU? + en da 0 6) 


imple two-dimensional, double-wedge airfoil. . . ae : = ‘ : ; : 
simple ¢ . 8 After differentiating Eq. (5) once and Eq. (6) twice with respect 


to /, substitution yields 


a@ + [la +e) + (co + s)U? — 20U/U)| a&@ + [ler + en 4+ eiez) + (exer + cies + 69 — C3¢9) U2? 4+ 
(2c, + c2)UU + cocyU! — (c, + 2c7)(U/U) + 2(U/U)? — (U U)la@ + [leper 4+ eyez) + (een + eyes) U2? + (eg 4+ Ges + 0007 
€x¢9) UU + B3eegUU3 — (ez + 2en)(U/U) + 2e,(U/U)? — (U/U) 4 cUUlea + 
lean + eycyl? + coc UU +¢6,UU — qe,(U/U -~ e,(U/U + 2e4,(U/U)2 Ila 0) (7 
The constants ¢, ..., ¢), are given below. 
c. = [4bpc/m|; co [(y + L)br2p/em|;) cz = [—(y + 1)rb2p/m); cs Wr, Cs C23 
(= [453 pc Sle); ts = [(4 T 1 )7*b3p Sl ac); cy = [ al de l th2p ee Cw C9; Ow Wa” 


Thus, the equations of motion have been reduced to a single fourth-order equation in torsion with coefficients which are functions of 
the forward flight velocity Ll’ 

Now, for convenience, assume L’ to vary linearly with time in a manner prescribed by the equation L’ at, where a is an acceleration 
of 10g. Also consider the airfoil to have the following parameter values: 


b = 1.5ft;7 = .04;m = 5.24slug/ft; 

Ie = 2.95 slug-ft?/ft; we 111.7 rad/sec; 
w, = 55.85 rad/sec; p = 0.00233 slug/ft3; 
Cc = 1,117 ft/sec. 


Using the above values, Eq. (7) becomes 


«a + [7.388 + 5.894 X 10>8t? — (2/t)] @ + [8.359 * 10-*%! — 17.6772 + 9.412 & 107% + 1.561 & 104 — 
(11.797/t) + (2/t2))@ + [2.508 & 10773 + 4.06222 — 35.355t + 5.092 * 104 — (2.497 & 104/t) + (8.818/t2)]a@ 4 
[—5.514 X 104%? + 2.965t + 3.892 * 107 — (3.717 X 104/t) + (2.495 X 104/t2)|a = 0 (8 


The coefficients of Eq. (8) were tabulated for various times ¢ 
near the expected critical point, thus yielding several fourth- 


order equations with constant coefficients. Routh’s stability re : a ee : , — 
criterion was then applied to the coefficients of each of the re- pect seer eqn ate ray Ss oe bore wapenree 5 — 
sulting characteristic equations, and the point at which a root with the critical point and apply ow Routh’s criterion to each of the 
a positive real part appeared was found to be at ¢ 17.86 sec resulting characteristic equations as before, the critical speed 
« . « < «< « a « = ‘a , . ° . . ~- . - c 
which corresponds to a speed of 5,750 ft per see as given by classical flutter analysis is found to be 5,645 ft per 
. ae ade . ‘ ° ° e - 2 
5 . —- : -c, Which gives a conservative e 0: ‘r sec, or 183 
Now, assuming negligibly small acceleration through the ores vhich eee + Sere error of } + Oo me, « , 
critical range (U U = 0), as is the case in classical flutter percent of the true flutter speed of 5,750 ft per nee? 
analysis, Eq. (7) with the same constant values as before becomes lt appears that the error introduced by neglecting the effects 
of acceleration is negligibly small and on the conservative side 


a@ + (7.388 + 5.894 X 10-"U2) a@ + and that classical flutter analysis may safely be used in investi 
(8.359 XK 107!2U4 — 17.677 XK 10°72U2 + gating the stability of systems moving through the critical range 
1.561 X 10%)a + (4.062 K 107~2U? + 5.092 XK 104)a@ + with a rapidly varying forward flight velocity. However, before 


(—5.514 X 10?U? + 3.892 X 10%)a 0 (9) this belief could be definitely established, numerous other cases 
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Fic. 1. Coordinate system and notation for an oscillating two 
dimensional airfoil 


over a wide range of accelerations would have to be considered, 
On the other hand, substantiation is given by Ref. 2, which gives 
the results of flutter tests on unswept rectangular planform wings 
propelled by ground-launched rockets, where the acceleration 
effects on flutter speed and frequency were found to be negligible 


for accelerations in the range 19-31g. 


REFERENCES 
Ashley, H., and Zartarian, G., Piston Theory--A New Aerodynamic Too 
or the Aeroelastician, Journal of the Aeronautical Sciences, Vol. 243, No. 12 
Dec. 1956 pp. 1, 109-1,119 
Molyneux, W. G., Ruddlesden, F 
Tests Using Ground-Launched Rockets, with Results for Unswep 


ARC R&M No. 2,944, 1956 
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On the Error in Plug- Type Calorimeters Caused 
by Surface-Temperature Mismatcht{ 


John C. Westkaemper 

Research Engineer, Defense Research Laboratory, and 
Assistant Professor of Aero-Space Engineering, 

The University of Texas, Austin, Tex 


y bea USE of plug-type calorimeters for measuring convective 
heating is common practice, due in part to the apparent 
simplicity of the method. On closer examination, however, 
proper use of the method is found to be less simple than it appears 
The problem of insulating the plug is difficult to solve in some 
cases,'!? and even an excellent insulator can induce errors in 
measurement Another large 
magnitude, that of surface-temperature mismatch, is the subject 


potential source of errors of 
of this note 

Rubesin® has specifically considered the effect of temperature 
mismatch on calorimeter accuracy, and concluded that very 
large errors could result if the calorimeter surface temperature 
were not closely matched with that of the surface upstream of 
the calorimeter. Using the nomenclature of Ref. 3, as shown 
in Fig. 1, Rubesin’s expression for the local heat-transfer 
coefficient of turbulent, incompressible flow with a step change 
in surface temperature is 


k(x, L 
(tw, — to)/(twe — to) and s = (tu — t tw. — to 


h(x, O)fb + s[1 — (L/x 


where 6 


7 This work was sponsored in part by the U. S. Navy Bureau of Weapons 


under Contract NOrd-16498 





oe Fic. 1. Flat plate nomen- 
clature 
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his equation was used to determine the ratio of the average 
Ww L) to 
the heat-transfer coefficient which would prevail at the center 
at a constant temperature 


heat-transfer coefficient over the calorimeter face 


of the calorimeter if the plate were 
The resulting relation is 


h(w, L 


(eho) 


Rubesin evaluated F(L/W) and H(L/W 
From Eq. (2) one may easily 


numerically and 


f the results 


included plots « 
calculate the ratio of the measured heat-transfer coefficient to 
the true value which would exist if there were no temperature 
mismatch 

In a later report,’ Reynolds, Kays, and Kline found that 
while Eq. (1) is valid for the region of test data available to 
Rubesin, a relation which holds over a wider range of Reynolds 
numbers is 


hix, L h(x, O){b + {1 L/x 9/0 ' 3 


The average heat-transfer coefficient to the calorimeter is 


oH +1) 
where O | q(x)dx J h(x, Lot toid 
Ji L 
So 
: | eu L \ 90 ’ 
hw, L J h(x,0)96 +2] 1 ( ) Sd 
W-—L I \ { 


which may be integrated to give 


hiw, L h(W, 0) 


j5 (1 — (L/W)*3)b 5 (L/W )(")" my ad, 
4 (1 —L/W $(1 — L/W)L\L { 


Numerica! values of 


H'( L/W 


Fic. 2 





\ modified form of Eq. (2) is then found to be 


h(w, L 


(= tL ) 
h 0 


where F(L/W 
modified term 
deviates from one by less than one percent The values of 
F(L/W) plotted in Ref 
Thus for L/W > 1/2, Eq. (4) becomes 


F(L/W 


is the same as in Eq. (2), but H’(L/W) is a 


W of 1/2 or more, F(L/W 


For values of L 


3 are erroneously high by a maximum 


of 12 percent 


h(w, L , f 
7 1+ H’(L/iU (ty, -f (t lo)} 5 
(" +L ) 
h ,0 
; ; 5 (L/W) W \ 9/1 8/9 
where H’(L/W) = : l l 
4 (1 L/W i 


For convenience, H’(L/W) is plotted as a function of L/W in 
Fig. 2; from this plot and Eq. (5), the performance of a calorim 
eter may readily be determined for any set of test conditions, 


subject to the restrictions inherent in Eq. (5 
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A comparison of Eq. (5) with the results of Ref. 3 indicates 
that the error predicted in Ref. 3 is excessively large. Rubesin’s 
original conclusion is still valid, however, as the revised results 
presented herein still indicate that temperature mismatch can 
result in very large errors in calorimeter results 


REFERENCES 
'! Westkaemper, J.C., An Analysts of Slug-Type Calorimeters for Measuring 
Heat Transfer from Exhaust Gases, AEDC TN 60-202, Nov. 1960 
2 Westkaemper, J. C., An Experimental Evaluation of the Insulated 
Technique of Measuring Heat Transfer at High Velocities, Rep. No. DRL 
The University of Texas, Austin, Tex., 
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139 
CF-2765, Detense Research Lab., 
Jan. 1959 

3 Rubesin, M. W., The Effect of an Arbitrary Surface-Temperature Varta 
tion Along a Flat Plate on the Convective Heat Transfer in an Incompressible 
Turbulent Boundary Layer, NACA TN 2345, April 1951 

4‘ Reynolds, W. C., Kays, W. M., and Kline, S. J., Heat Transfer in th 
Turbulent Incompressible Boundary Layer. II—Step Weall-Temperature 
Distribution, NASA Memo. 12-2-58W, Dec. 1958 


Magnetogasdynamic Effects on Laminar 
Flame Propagation 


James V. Beck* and Tau-Yi Toong** 
Avco Research and Advanced Development Division, 
Wilmington, Mass. and M.1.T., Cambridge, Mass., respectively 


March 14, 1961 


A THEORETICAL STUDY of the propagation of a steady, one- 
dimensional, adiabatic, laminar flame in a chemically reac 
tive gaseous mixture in the presence of a uniform, transverse, 
magnetic field has been conducted at the Massachusetts Insti- 
tute of Technology. This note briefly reports some of the re- 
sults obtained in this investigation 

The chemically reactive gaseous mixture under consideration 
is assumed to be an electrically conducting (though electrically 
neutral) continuum within the laminar flame. For an ob 
server standing on the steadily propagating flame, the con- 
servation equations of mass, momentum, and energy, after com- 


bination with Maxwell’s equations, are as follows: 


(d/dx}(pu) = O (1 
du d B 1d du 
pu + 2 =| —_ pv = 0 (2) 
dx dx 2u 3 dx dx 


aT d ( «Be? 4 du\* 
ata dx as dx ii Qu 3 dx : 


d dT d {B 
r + u,B,, —q 0 (3 
dx dx dx \u 


Furthermore, one obtains the following relationship from Max- 
well’s equations and Ohin’s law: 


(d/dx)(B/p) = o (uB — u,B,) (4) 


In the above equations, p, p, T, cp, v, 4, w, and o are the mass 
density, mechanical pressure, absolute temperature, specific heat 
at constant pressure, kinematic viscosity, thermal conductivity, 
magnetic permeability, and electrical conductivity of the gaseous 
mixture, respectively, « is the mixture velocity in the direction 
of the flame propagation (x-direction), B is the magnetic-field 
strength in the direction of the imposed magnetic field (normal 
to the x-direction), g is the energy source due to chemical reaction 
per unit volume and unit time, and the subscript , designates 
the conditions upstream of the flame. It is to be noted that, 
in the derivation of the energy equation, the specific heats of 
the various components of the mixture are assumed to be equal 
to one another. 

The boundary conditions are the following: For x = x,: 


* Staff Scientist 
** Associate Protessor of Mechanical Engineering. 
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Fic. 1. Magnetogasdynamic effect on burning velocity. 


p Pur P = Pw B B 


(5 
du/dx = dT/dx 0 
For « = x, (at the completely reacted condition 
du/dx = dT/dx = 0 (6) 


Examination of the above differential equations and their 
boundary conditions indicates that this is an eigenvalue prob- 
lem—i.e., a solution can exist only for an eigenvalue of the flame- 
propagation speed u,, which is usually called the burning ve- 
locity and designated as S, 

An integral method is employed in this study to investigate 
the effects of a transverse magnetic field on the burning velocity. 
The usefulness of this method has been well demonstrated in 
boundary-layer problems, in problems where one is interested 
only in overall effects. In laminar-flame-propagation problems, 
Spalding! has further shown that, with a linear temperature 
profile in the flame zone, the approximate solution obtained by 
the integral method agrees within 20 percent of the exact solu- 
tion. It is expected that the accuracy of this method would not 
be greatly impaired by the added complexities of the present 
problem due to the presence of electromagnetic action 

Fig. 1 shows the effect of a transverse magnetic field on the 
steady propagation speed of an adiabatic laminar flame, for the 
case where the reactive mixture is assumed to be a perfect gas 
with constant specific heats, thermal conductivity, magnetic 
permeability, and electrical conductivity, and where viscous 
dissipation and kinetic energy are negligible. Therein, the 
ordinate ¢ is a dimensionless quantity, containing the burning 
velocity S,, thermal diffusivity \/p.c,, temperature rise (7) 


T,,), thermal conductivity \, and the integrated volumetric rate 
*Td 


of chemical-energy change with respect to temperature gdT; 
Tu 
the abscissa 8 is the ratio of the magnetic enthalpy B,?/up, 
to the sensible enthalpy c,7\,, both corresponding to the un- 
reacted condition; and the parameter a is the magnetic Prandtl 
number, which is defined as the ratio of the magnetic diffusivity 
(uo)! to the thermal diffusivity A/pucp. The curves shown are 
for specific values of the ratio of the specific heats y and of the 
temperature ratio 7/7, or (T4/T.)o [where (74/7T.)o is the tem- 
perature ratio when there is no imposed magnetic field], and are 
obtained by the use of a linear profile for the temperature dis- 
tribution and a hyperbolic profile for the distribution of the 
magnetic-field strength within the flame. When, however, o 
is infinite, it is no longer necessary to assume this latter profile, 
because it can be obtained from the familiar relationship that 
the product “1B remains constant throughout the flame. For 
such a case, the B-profile can be approximated very well by a 
hyperbola. This finding leads to the choice of a hyperbolic 
B-profile for the general case of finite electrical conductivity 
Several observations can be made from Fig. 1. First, the 





— xe hlUrlUrM 


> ~~, OS 








(6) 


their 
prob- 
flame- 
ig ve- 


tigate 
ocity. 
ed in 
rested 
lems, 
ature 
-d by 
solu- 
d not 


esent 


1 the 
r the 
t gas 
netic 
cous 

the 
ning 


rate 
dT; 


Pu 
un- 
ndtl 
vity 
are 
the 
em- 
are 
dis- 
the 
» Co 
file, 
hat 
For 
ya 


lic 


the 





READERS’ 


effect of the magnetic field is to increase the nondimensional 
burning-velocity parameter ¢. Secondly, there is an asymptotic 
value for this parameter as the enthalpy ratio 8 approaches in- 
finity, this asymptotic value being equal to the ratio of the spe- 
cific heats, regardless of the value of the magnetic Prandtl num- 
ber a. Thirdly, for large values of a, g attains a maximum as 
8 increases. Fourthly, contrary to the expectation that an in- 
crease in the electrical conductivity would result in an increase 
in the magnetogasdynamic effect, the results of this analysis 
indicate that the burning-velocity parameter decreases with 
increasing value of ¢ at a given enthalpy ratio.* 

The temperature rise (7; — 7.) across the laminar flame can 
be computed by integrating Eq. (8) between x, and x,. For a 
given change in the chemical composition, it can be shown that 
this temperature rise is independent of the magnetic Prandtl 
number, increases with increasing enthalpy ratio 8 and approaches 
an asymptotic value as 6 approaches infinity. In particular, 
for the case under consideration, the maximum ratio of the tem- 
perature rise with and without magnetic effect is equal to the 
ratio of the specific heats. Thus, if 7, is maintained constant, 
the maximum value of 7),/7,, is 15 for y = 1.4 and (7;/7,)o = 11. 

The two limiting curves for 7,/7, = 11 and 15 are shown in 
Fig. 1 for a specific value of the magnetic Prandtl number, to- 
gether with the curve for (7),/7,,)y) = 11 and 7, remaining con- 
stant as the enthalpy ratio varies. As expected, this last curve 
approaches the two limiting curves at low and high values of 
the enthalpy ratio 

Notwithstanding the increase in the temperature rise in the 
presence of a transverse magnetic field, it is expected that the 

1b 
much greater increase in the value of gdT with increasing 
Tu 
temperature, coupled with an increase in the value of ¢, would 
result in an increase in the burning velocity for a given chemically 


reactive mixture 


REFERENCE 
1 Spalding, D. B., Approximate Solutions of Transient and Two-Dimen 
sional Flame Phenomena: Constant-Enthalpy Flames, Proc. Royal Society, 


A-245, pp. 352-372, 1958 


* It is to be noted that the curve for an infinite value of @ represents an 


upper bound for the magnetic effect as o> 0. It is different from that for 


which ¢@ is identically zero. For the latter case, it is obvious there is no 


magnetic effect at all 


Some Detrimental Aspects of the Evaporation 
of a Quartz Shield for a Particular Re-Entry 


Problem 


Ernst W. Adams 

Aeroballistics Division, G. Marshall Space Flight Center, NASA 
Huntsville, Ala. 

March 22, 1961 


c be TRANSIENT PERFORMANCE of a foamed-quartz shield 
is treated herein for the spherical and the conical surfaces 
of a vehicle (Fig. 1) which weighs 8,640 kg, has a ballistic 
factor of 500 Ib/ft?, and re-enters on a trajectory defined by the 
altitude H above sea level and the speed V as functions of the time 


t (Fig. 1). The re-entry angle at H = 120 km is 94.5°; a 
maximum deceleration of 7.7 g is reached at H = 24 km. The 


following properties are assumed for the quartz shield: a thermal 
conductivity of 2.5 X 1074 keal/(m°K sec) a specific heat of 


0.29 keal/(kg°K) a specific weight y = 1,050 kg/m‘, a heat of 
vaporization h, = 3,050 kcal/kg, the measured functions yu(7) 


and p,(7) for the viscosity and the vapor pressure, and a constant 
temperature of T = 300°K at H = 120 km. 

The mathematical solution for the heating and the melting of 
the shield includes the temperature 7 = T7{(x,y,t) and the ve- 
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Fic. 1. The re-entry vehicle and its trajectory 
locity components u = u (x,y,t) and v = v(x,y,t) as functions of the 


time ¢ and of the surface-fixed coordinates x and y. Refs. 2-5 
are employed for the aerodynamic heating of the shield and Ref 
1 for the evaporation of the molten quartz. Even at a tem- 
perature of 4,000°K, quartz is so viscous that the inertia forces 
of a quartz flow, caused by the external forces, are negligible 
For a small vicinity of the cross section x = 0 at the stagnation 
point of the vehicle, Ref. 6 presents a numerical method which 
vields exact transient solutions for 7, u/x, and v. One unified 
system of differential equations is valid for the total shield, 
x > O, if the dependencies on x of 7, u/x, and v are neglected 
for x > 0. The solutions show for the total shield that (1 
OT/dx| < 0.015 |OT/dy|, (2) |dv/Ox| < 0.02 |dv/dy|, and 
(3) there is practically no flow of molten quartz because it 
evaporates so readily. Thus, the derived solutions are good 
approximations to solutions which do not rest on the said con 
ditions. 

The necessary thickness d*(x) of the shield is determined so 
that 7(x,d*,t) approaches an upper limit of 460°K as 07(x,d*,t) + 
ov > 0 at the supposedly insulated inner edge of the shield. The 
table presents d*(x), the time-integrated ablation s(x), the initial 
weight yd*(x), the ablated weight ys(x), the evaporated weight 
W(x), the total aerodynamic heating rate Qsero(x), the corre- 
sponding figure Qsero(x) for suppressed vaporization, the total 
amount of heat Qpi(x) = Quero(®) — Qaero(x) + h, W(x) blocked 
by the evaporation, and the total amount of heat Qraa(x) which 
is radiated into the air from the surface of the supposedly opaque 
shield, the emissivity constant of which is assumed to be « = 0.4 
The absorption and the emission of radiation in the air boundary 
layer are neglected 

The table shows that Qpi(0)/Qaeo(O) = 0.71—Le., 71 percent 
of Queso is blocked by the quartz vapor from entering the shield at 
x = 0. The term Qsero(O) — Qaero(O), which represents the effects 
of the vapor diffusion across the air boundary layer, accounts for 
65 percent of this blocking action. At the station » 4(m), 
Qrad(4)/Qaero(4) = 0.80, which means that 80 percent of Qaero 15 
radiated back into the air. On the spherical nose, 0.98 < 
W(x)/ys(x) < 1.0, or the ablation consists practically of evapora 
tion only. The table’s bottom line gives the integrals of the data 
presented per unit area over the spherical and conical surfaces 
of the vehicle. The shield loses only 2.7 percent of its initial 
weight by ablation. This blocks, however, 16.9 percent of the 
aerodynamic heating rate S sQuero(% 2ar(x) dx, 71 percent of 
which is radiated back into the air for the entire shield 

Detailed investigations show that the evaporation of the 
quartz shield is very sensitive to increases in Quero(X) When the 
heat transfer coefficient h(0,t)kg/(m? sec) at x = O is increased 
by a constant factor of 2, s(0) rises by a factor of 2.43, and, due 
to a reduced thermal penetration of the shield, d*(Q) rises by a 
factor of 1.65. When h(4,t) is increased by this factor of 2 at 
x = 4m, where s(4) is negligible, d*(4) rises only by a factor 
of 1.05. In addition to these detrimental consequences of the 
evaporation of the quartz, it is seen in the table that d*(Q) at 
0 can be reduced by a factor of 0.84 when the evaporation 
The nonevaporating shield, which 


x= 
of the quartz is suppressed 
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TABLE | 
Vaporizing Quartz Nonvaporizing Quartz 
\ a” ys yd* W Du So Oy Ora , d* iors QO \ 
~ % SI 
m1 mim mm kg/m? kg/m? kg/m? keal/m? keal/m?  keal/m? keal/m? m mm mm keal/m?  keal/m? ; 
0 38.45 75.90 40.37 79.69 40.23 $92,089 263.402 351,395 135,037 O 16.57 63.838 472.015 446.486 ‘ 
0.25 30.15 69.85 31.66 73.34 31.46 107.201 227.579 275.566 125.472 a 
0.50 9.90 58.61 10.39 61.54 10.26 192,619 131,751 92,147 92 076 
1.50 1.30 54.52 1.36 5/25 1.36 80,605 72.103 12.646 59,167 re 
2.50 0.74 353.938 0.78 56.63 0.77 70,509 66 O89 6,768 54,902 
$00 0.45 53.25 0.48 55.92 0.47 61.887 60.007 3,311 19 833 
Surface 
Integrals 17.04 1707 16.60 2,360,190 2,102,661 399.662 1,696,678 
kg kg kg keal keal keal keal 
loses s(0 16.6 mm. of its initial thickness at a due to the flat-earth system. Therefore, the component of longitude rat¢ ; 
flow of molten quartz, reaches a maximum surface temperature in the radial direction, @ sin ¢, was neglected. The inclusion C 
of 4,623°K and radiates 94.6 percent of Qner(O) back into the of this term does not unnecessarily complicate the problem 
air. The evaporating shield reaches only 3,101 °K and radiates By rearranging Eq. (3) it is seen that 
27.5 percent back. It may be concluded for the re-entry problem M 
at hand that the prevention of the evaporation of the quartz dy ( a sin adv) mg/CyS)VV| — sin dé j R 
would allow a thinner quartz shield which could cope with even ~ * 
considerably increased heating rates Eq. (4) can be further simplified by letting - 
V2/¢ n a r 
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(9) now becomes | 
i cE do L, D)(cos B/2) sin ~ldyn/(1 n 10 st 
— Eqs. (8) and (10) cannot be solved analytically for the constant F 
An Improved Method for Determining the ' oe , 7 , ° rene ti 
Papi - yank-angle case owever, very good results can be obtainec 
Lateral Range of a Gliding Entry Vehicle , 
by holding latitude and deflection angle constant and integrating ? 
over small increments of 7. The following approximations result 
W. Scott Jackson l 
Aerodynamicist, Douglas Aircraft Co., Inc., ; ; if L in al In MN a SA 
F ? 1 4 Ss i “OS f é py xX 
Missiles and Space Systems Engineering, en wN NAD me am PN ‘ 
Santa Monica, Calif , 
March 24, 1961 — | nN 4 ( 
cos wyf In 11 
1 —n, i 
A RECENT publication by Slye! indicates that the analytical 
estimates of lateral range (or latitude, ¢, in the case of 
equatorial orbits) of gliding entry vehicles tend to be very a 
optimistic at large ranges. The reason for this is evident from | 
the simplified equilibrium equations of motion for a vehicle at “i g 
very small flight-path angles and a constant bank angle, 8 = 
zm 4; ~ 
| ° FLAT EARTH 4f 
rosy ro (METHOD OF 
-~ pl?C,S cos B meg + (mV2/) 0) (1 8 REFERENCE |) . , * 
a a. OH ; r Fic. 1 Comparison of lat me 
‘ l : = y eral range estimates 
m4 -— pV2CpS (2 & | ’ _ 
» 2} & 
he a 
4 rs | t] 
= ec l cad Tee ~ | L SPHERICAL EARTH 
mViy + @sin @ = pV2C_S sin B (3 < | | J7 tt 
S Oe 
i in 
Eqs. (1) and (2) are straightforward; however, Slve assumed ol. a di 
. : ts) ' 2 3 4 
that the velocity vector was deflected through an angle, y, in a LIFT-DRAG RATIO la 
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\ comparison of the results of this method with the method of 
Slve is shown in Fig. 1 for a bank angle of 45° and an » at A 0) 
of 0.9. The integrations are terminated when the velocity 


vector is pointed due north—i.e., when the deflection angle is 


YQ) \s expected, the change from a flat- to a spherical-earth 
ssumption in Eq. (3) causes a considerable difference in the 


results at large ranges 
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Improved Numerical Solution of the Blasius 
Problem With Three-Point Boundary 
Conditions 


William J. Christian 

Research Engineer, Armour Research Foundation or 
linois Institute of Technology, Chicago, III 

March 28, 1961 


= VELOCITY distribution resulting from laminar, constant 
pressure mixing of a fluid laver at rest and a moving stream 
of the same fluid in steady two-dimensional flow has been de 
boundary 


scribed by the Blasius equation with three-point 


conditions, that is 


f ff 0) ] 
(0 0) Ya 
lim {’(y () (2b 
n— 
lim f’(y 2 2c 
_—> 


where n and f(y) are the usual Blasius variables 

Solutions have been given for this problem,? us well as 
for more general cases where the two layers consist of different 
fluids or where both are in motion © In addition to their 
importance in the familiar jet-mixing problem, such solutions 
ilso have application in an approximate way to the analysis of 
separated and wake flows,’:> and accurate numerical results 
should be of interest This note is concerned with the par 
ticular problem defined by Eqs. (1) and (2) and the achievement 
of higher accuracy in its numerical solution 

In connection with a numerical procedure programed for th« 
UNIVAC 1104, a routine was developed for generation of the 


) 


solution of Eq. (1) with boundary conditions given by Eqs. (2 


{ 


for arbitrary increments in 7. The procedure is based on analytic 
continuation of the function through Taylor’s series, and pre 
supposes that the initial values f(0), (0), and f"(O) are known 
For this reason, the method is not of direct value in the numerical 
a three-point boundary-value problem involving a 


With only 


minor modification, however, the computer program provides 


solution of 
nonlinear differential equation such as Eq l 
i means to improve the accuracy of previously calculated values 
ff'(O) and f”’(0), and thereby to generate the complete numerical 
olution with six-decimal accuracy 

The concept of analytic continuation is not new, of course, 
but it has been most useful as an analytical tool rather than as a 
\vailability of 


computation now makes it feasible to use the method for numeri 


numerical procedure high-speed automatic 


cal integrations of single-point boundary-value problems such 
that, within the limits of convergence of Taylor’s expansion, 
truncation error may be made arbitrarily small regardless of the 
interval size. Thus, numerical integrations whose accuracy 
depends only on round-off errors may be developed for rather 


large increments in the independent variable. A brief descrip 


FORUM $11 








TABLE | 
The Function f(7) and its Derivatives for —14 Z » Z 4 
7 f(m) {(m) ™(n) 
-14.0 -1.238495 0.000000 0.000000 
-12.0 -1. 238494 . 000001 . 000001 
-10.0 -1.238488 . 000008 . 000010 
- 8.0 ~1.238413 . 000101 .000125 
- 7.0 -1.238214 . 000347 . 000430 
--6.0 -1.237527 . 001198 . 001 484 
- 5.0 -1.235157 . 004131 -005109 
- 4.4 -1.231483 . 008672 . 010709 
- 34 -1.223776 . 018175 .022376 
- 3.2 -1. 207650 . 037963 . 046433 
- 2.8 -1.188075 . 061810 . 074998 
- 2.4 -1. 156284 . 100134 - 119922 
- 2.0 -1. 104988 - 160918 . 188653 
- 1.8 -1. 068740 . 203092 - 234492 
- 1.6 -1.023082 . 255301 . 289101 
- 1.4 -0. 965831 - 319340 . 352793 
1.2 -0. 894438 . 396991 . 425033 
- 1.0 -0. 806022 . 489804 . 503973 
- 0.8 -0.697437 . 598787 . 585948 
- 0.6 -0. 565422 . 724000 . 665098 
- 0.4 -0. 406840 . 864104 . 733354 
- 0.2 -0. 218992 1.015975 - 781112 
0 0 1.174540 . 798828 
0.2 . 250821 1. 333022 . 779453 
0.4 . 532688 1. 483706 . 721076 
0.6 . 843284 1.619168 . 628665 
0.8 1.178951 1. 733663 . 513761 
1.0 1.535145 1. 824200 . 391760 
1.2 1.907037 1. 890909 . 277731 
1.4 2.290102 1.936564 . 182562 
1.6 2. 680549 1.965524 . 111066 
1.8 3.075515 1.982525 . 062463 
2.0 3, 473043 1.991751 . 032452 
2.4 4.271431 1.998523 . 006896 
2.8 5.071171 1.999304 . 001064 
3.2 5.871140 1.999981 .000119 
3.6 6.671137 1.999999 .000010 
4.0 7.471137 2.000000 . 000001 





tion of the application of analytic continuation to Eq. (1) is 
given below 

In order to effect analytic continuation, the function f(»), 
‘(n » and f” yn) are represented by a Taylor’s series expansion 


about any ny as follows 


where {"(n SS” nln lia, (7 ) ” ) 
rw 
) 

‘ n n n! 6 


Now, if the differential equation |E 1)| is to be satisfied, a 


recurrence relation exists among the coefficients a, such that 


n—1 


where 


By means of Eqs. (6) and (7), it is possible to obtain the 


Taylor’s series |Eqs. (3), (4), and (5)| to any desired number of 


terms at any m for which f, f’, and f” are known Accordingly, 
f(no + An), f’ (no + An), f'(n + An), f"(n + An) can be found 
for an interval Ay within the radius of convergence of the series 
With the exception of round-off error in the computer, any 
desired accuracy may be achieved by specifying the maximum 
truncation error in the Taylor expansion. The maximum value 
which Ayn may have is determined by the radius of convergence 
and by the maximum rounding error that can be tolerated. If 


An is too large, individual terms of the series (which contains 
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both positive and negative terms) may be many orders of magni- 
tude larger than the sum of the series, and rounding errors will 
be magnified. 

The procedure for generating values of the function and of its 
derivatives is thus straightforward, as is any numerical integra- 
tion starting with single-point boundary conditions. Applica- 
tion of the procedure to improvement of the accuracy of f’(0) 
and f"(0) involved a simple iteration technique, starting with 
published values of f’(0) and f"(0). This consisted of a number 
of machine runs for various values of f’(0) and f’(Q) to provide 
numerical relationships between f’(—  ), f’(@), f’(0), and f”(0) 
Parabolic interpolations were then used to establish f’(0) and 
f"(O) such that, within eight significant figures, 

lim f’(n) = O 


> — « 


lim f’(n) 
o> 0 
Although computations were carried out in a manner designed 
to produce eight significant figures in the results, round-off errors 
in the UNIVAC 1105 make the eighth figure (and possibly the 
seventh figure) uncertain. To check the effect of round-off 
error, calculations were performed using increments An of 0.1, 


Inertial Guidance Systems 


SCIENCES NOVEMBER 1961 


0.2, 0.5, and 1.0. Comparison of the results indicates that 
accuracy to at least six decimal places has been achieved in the 
function f(y) as well as in its first and second derivatives, 
Numerical results are given in Table 1. 
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(Az), — H(¢d: + Qe sind + ¢$,2,) 


which was used in relations (1). 


Appendix B 


Sample Calculation of Alignment Time 


Second-Order System 

Assume that a; = 0.lg at | cps and ¢,(0) and ¢.(0) = 
100 mr. The accuracy required is ¢, = 10 mr. 

(1) From Fig. 7,7 = 170sec. 

(2) From Fig. 8, curve B,, 3.8 r’s required to reduce 
¢: by a factor of 10. The time for closed-loop align- 
ment is 660 sec. 

(3) From Eq. 8, (¢,)pre = 170/(2 & 10‘) (100 mr) = 
0.85 mr. 

(4) From Eq. 9, 7) = [0.1 g/(0.85 X10~*) rad] X 1 


27 = 19sec. 


(5) From Step (3), ¢, must be reduced by a factor 
of 100/0.85 = 120 and from Fig. 8, curve A, this will 
take 57's or from step (4) 95 sec. 

(6) The total alignment time is 660 + 95 = 755s 
or 13 min. 


Appendix C 


Inherent Filtering of Velocity Meter 


Because the output wheel of a velocity meter turns 
at a rate proportional to acceleration, the counter on 
the output wheel indicates the integral of acceleration 
(Sketch A). 

This instrument can thus be converted, via a simple 
feedback loop (to its torquer), to act as an acceleration 
filter of any desired time constant (Sketch B), which 
converts the block labeled ‘‘Accel.”’ in Fig. 4(b) to the 
one labeled ‘Accel. with filter” in Fig. 10. 








Microcards for Volumes 1-5 Available 


Microcards for Volumes 3-5 of the JouRNAL are now available (Vols. 1 and 2 were republished in this form 
earlier). All five volumes are priced at $5.50 each, and are obtainable through: 


J. S. Canner & Company, Inc. 
Boston, Roxbury 20, Mass. 
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makes it difficult to establish the identity of the author. This is 
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name of the organization with which the author is associated 
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which the paper is received will be inserted by the Editor. The 
author’s title or position should be indicated in a footnote 

SUMMARIES OR ABSTRACTS: An abstract to be printed at the 
beginning should accompany each article. It should be adequate 
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major conclusions reached, since summaries in many cases con- 
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tific reference indexes. Abstracts printed in other journals, espe- 
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the major conclusions to a nonspecialist in the subject and should 
contain from 100 to 300 words, depending on the length of the 
paper. 

SUBHEADINGS: Subheadings should be inserted by the author 
at frequent intervals. The work of editorial preparation will be 
simplified by the author’s provision of many subheadings. 
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materials or processes or of preliminary experiments or theories 
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REFERENCES: References should be numbered consecutively 
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the part to be reproduced. Each table should have a caption 
such as Table 1, Table 2, Table 3, etc. Captions should be com- 
plete in themselves so as to make the data intelligible to the 
reader without reference to the text. A duplicate list of captions 
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Use “‘Fig. 1’’ (not Figure 1), ‘‘Figs. 3 and 4,” etc., in both the text 
and the numbering of illustrations. In the text, “Eq. (1)” or 
“Eqs. (1) and (2)” is used, not “Equation (1).”" In captions, 
legends, and in table headings, write all words in full; do not 
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MATHEMATICAL WoRK: Formulas may be typewritten or 
carefully written in pen and ink, the writing to be large enough so 
that it can be marked for the printer. Considerable space for 
marking should be allowed above and below all equations. All 
complicated equations should be repeated on separate sheets with 
adequate space left for marking. The solidus should be used for 
simple fractions appearing within the text. Make all expressions 
ciear to the typesetter. Greek letters used in formulas should 
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SPERRY GYROSCOPE COMPANY, DIVISION OF SPERRY® 
RAND CORPORATION 


STANDARD-THOMSON CORPORATION 
STANLEY AVIATION CORPORATION 
THERM, INCORPORATED 

THIOKOL CHEMICAL CORPORATION 
THOMPSON RAMO WOOLDRIDGE INC. 
TRANS WORLD AIRLINES, INC. 


UNION CARBIDE METALS COMPANY, DIVISION OF® 
UNION CARBIDE CORPORATION ‘ 


UNITED AIR LINES, INC. 

UNITED AIRCRAFT CORPORATION 

UNITED STATES AVIATION UNDERWRITERS, INCH 
VEHICLE RESEARCH CORPORATION 

WESTERN GEAR CORPORATION 


WESTINGHOUSE ELECTRIC CORPORATION 




















